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I .  INTRODUCTION 


The  objective  of  the  present  investigation  is  to  develop  a  finite  element 
model  and  computer  program  for  the  purpose  of  handling  elastic-plastic 
material  with  rate-dependent  yield  conditions.1  The  effort  is  in  two 
parts.  The  first  part  involves  the  analytical  development  in  which  the 
appropriate  incremental,  stress-strain  relations  are  developed.  The  second 
part  of  the  effort  involves  developing  a  finite  element  computer  program 
which  incorporates  the  analytical  development.  This  computer  program  will 
be  based  on  previously  developed,  approximate  three-dimensional  elastic- 
plastic  computer  code,  SANX.2’3  The  code  SANX  is  designed  to  perform 
structural  analysis  on  cylindrical  configurations  which  are  approximately 
axisymmetric  and  which  have  definite  nonaxi symmetric  features.  The  code 


SANX  was  developed  both  for  elastic  and  elastic-plastic  materials  with  no 
time  dependent  properties.2*3 


The  starting  point  for  the  development  of  a  three-dimensional  finite 
element  code  is  the  one-dimensional  analysis1  which  formulates  the  visco¬ 
plastic  response  in  terms  of  the  effective  plastic  strain.  This  approach 
is  extended  to  the  development  of  the  three-dimensional  model  in  the 
present  investigation. 


1.  W.H.  Drysdale,  "A  Theory  of  Rate  Dependent  Plasticity,"  Ballistic 
Research  Laboratory  Report,  APG,  MD.  (Forthcoming) 

2.  A.R.  Zak,  J.N.  Craddock  and  W.H.  Drysdale,  "An  Elastic-Plastic  Analysis 
of  Non-Axisymmetric  Structures,"  International  Journal  of  Computers 
and  Structures,  vol.  10,  pp.  841-846,  1979. 

3.  J.N.  Craddock  and  A.R.  Zak,  "An  Approximate  Finite  Element  Method  of 
Stress  Analysis  of  Non-Axisymmetric  Bodies  with  Elastic-Plastic 
Materials,"  Technical  Rept.  UILU-ENG  79  0501,  Aeronautical  6  Astronautical 
Engineering  Dept.,  University  of  Illinois,  Urbana,  March  1979. 
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II.  RATE  DEPENDENT  MATERIAL  MODEL 


The  rate  dependent  plasticity  model  was  introduced  for  isotropic 
material  and  used  in  the  analysis  of  uniaxial  stress  case.1  The  same 
model  will  be  used  in  the  present  investigation  with  a  slight  modification 
to  allow  for  use  with  orthotropic  materials^' when  needed.  The  yield  _ 
condition  for  orthotropic  materials  will  be  represented  by  Hill's  criterion^ 
and  this  reduces  to  the  octahedral  shear  stress  criterion  in  the  limit  for 
isotropic  materials. 1 

Using  the  rate  dependent  model  for  the  yield  criterion*  the  yield 
function  is  taken  in  the  form: 


f(o. .  -  a.  .)  =  K(e? .) 


(1) 


where  a.,  represents  the  strain  hardening  parameters.  The  rate  dependence 
is  defined  by  the  function  K  which  is  represented  by: 

K(e?j)  =  [1  ♦  blnfl  ♦  zP/zoH2  (2) 

Equation  (2)  is  an  empirical  formula  which  contains  material  constants  b 

and  e  .  The  dependence  of  K  on  the  rate  is  through  the  variable  ep  which 

will  be  defined  later  as  the  effective  plastic  strain  rate.  It  may  be  noted 
that  for  isotropic  material1  the  function  K  in  equation  (2)  is  multiplied  by 
the  square  of  the  uniaxial  yield  stress  at  zero  rate.  In  the  case  of  the 
orthotropic  material  there  are,  in  general,  six  yield  stresses  and  it  is 
not  possible  to  separate  one  stress  from  the  yield  condition.  The  six  yield  2 

stresses  are  included  in  the  function  f  on  the  left  hand  side  of  equation  (2) . 

Using  equation  (2)  to  represent  rate  dependent  yield  condition  for  orthotropic 
materials  implies  an  assumption  that  each  yield  stress  is  dependent  on  the 
rate  by  the  same  relation  to  the  effective  plastic  strain  rate. 


The  next  step  in  the  development  is  to  obtain  an  incremental  stress- 
strain  relation.  From  the  plastic  flow  rule  the  plastic  strain  changes 
are  related  to  derivative  of  yield  function2; 


(3) 


Using  the  definition  of  kinematic  strain  hardening: 

da..  =  Cde?.  =  Cdx||—  (4) 

xj  ij  Sc.j 

2 

where  C  is  the  strain  hardening  parameter.  During  the  incremental  plastic 
deformation  the  stress  and  strain  changes  must  remain  on  the  yield  surface 
and,  therefore,  from  equation  (2): 


M —  do. .  ♦|L-da1. 
9oij  1J  9aij  1J 


3K 

9eij 


de?. 

il 


(5) 


6 


(6) 


It  can  be  shown  that: 


3f 

*« 


3f 

3o7 


13 


and,  therefore,  combining  equations  (4)  and  (6)  with  (5)  gives: 


S^d0u  •  CliX  3377  5577  •  ^ 


ij 


'ij  ""  "ij  ’  ij 

Solving  for  the  parameter  dX  from  equation  (7) : 

1 


■*%  ■  0 


dX 


f3f 


«ii 


Combining  equations  (3)  and  (8)  and  changing  the  repeated  indices  in 
equation  (8),  for  clarity,  gives  the  incremental  change  in  the  plastic 
strain : 


(7) 


(8) 


dtf. 

13 


5f~ 


f3f 


3okl  d0nm  "  3eP 


— —  dep  1 
3eP  H 
mn 


3f 


(9) 


Consider  now  the  elastic  stress-strain  relation  for  incremental  changes: 
■  B^Cdet,  -  dePj) 


ij 


ijklv  kl 


(10) 


where  Eijkl  represents  elastic  material  properties  matrix  and  de^i  the 
total  strain  changes. 

Returning  to  equation  (7)  and  substituting  for  do. .  from  equation  (10)  and 
rearranging  gives :  J 


d>  ■  "tUrr E 


.de. 


3K 


ij  «kl  kl  3% 


de^) 


(11) 


where  0  is  by  definition: 


n  .  rr  3f  3f__  .  3f  -  3f__,-l 

D  30^  fcijkl  3oklJ 


(12) 
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Substituting  equation  (11)  into  equation  (10)  gives: 

do  •  •  s  [  E .  . .  -  -  DE .  .  E  ,  .  *  *■—  ■  ]  d£,  . 

xj  1  ljkl  ljrs  mnkl  30^  3 o^J  kj 


♦  DE. 


3f  3K 


ijrs  3o~  gjp”  °  kl 
rs  3ekl 

2 

Using  shorthand  notation  previously  introduced  for  elastic-plastic  materials 
permits  writing  equation  (13)  in  a  short  form: 

.  .  ,  __  3f _  3K  ,«p 

doij  “  Aijkldekl  +  DEijrs  3o  a*p  dekl  (14] 

rs 

It  may  be  noted  that  if  the  strain  rate  term  is  neglected  on  the  right  hand 
side  of  equation  (14) ,  then  the  remaining  terms  represent  the  relation  be¬ 
tween  total  incremental  stress  and  strain  changes  for  elastic-plastic  material 
used  in  SANX  model.  The  next  step  in  the  development  of  the  rate  dependent 
model  is  to  introduce  the  effective  plastic  strain  increment  defined  by: 

dep  -  ^2/3  deP.dePj  (IS) 

The  objective  of  this  will  be  to  use  this  concept  to  replace  the  strain  rate 
term  on  the  right  hand  side  of  equation  (14).  Using  flow  rule,  equation  (3), 
in  equation  (IS)  results  in: 


de*  =-  dX  *2/3 


^ij  ®°ij 


Returning  to  the  yield  function  and  the  function  K,  an  incremental  change  in 
this  parameter  can  now  be  written  as: 

_  3K  _  3K  3ep  ,ip  ,, 


Therefore: 


«_  .  2£_  diP 

SePj  l>  )cp 


uta'.  - 


? 


Returning  now  to  equation  (7)  and  substituting  for  dX  from  equation  (16)  and 
using  equation  (19)  results  in: 


3f  „  3f  3f 

3aij  «  '  30ij  30ij 


2/3 5577  — 

ij  iJ 


The  second  term  in  equation  (20)  can  be  simplified  and  equation  rearranged  as 
follows: 


15-  a?  .  3/2  c  '  2/3  H-  If-  d£ 

9?  Mij  a°ij 


P  3f 

=  aatj  d0ij 


Equation  (21)  can  be  compared  directly  to  equation  (III. 5a)  of  Reference  1, 
for  isotropic  materials,  with  the  following  substitutions: 


*L-2/3c 
3eP  y  3eP 


2/3  s  2/^y  ( 

"ij  "ij  y 

Before  proceeding  further  with  equation  (21)  it  is  useful  to  obtain  the 
expression  for  rate  of  change  of  K.  Using  equation  (2)  and  differentiating: 


2[1  ♦  bln(l  ♦  7~)1  re¬ 


consider  now  equation  (21)  and  apply  it  to  an  interval  of  loading  over 
which  the  rate  of  change  of  stress  is  represented  by  a  constant. 

General  equation  (21)  has  variable  coefficients.  However,  if  it  is  applied 
to  a  small  interval  of  loading  the  coefficients  can  be  assumed  to  be  constant 
over  this  interval.  Over  such  interval  the  stress  variation  can  be  approximated 
by  a  linear  variation  with  time: 

A^ij  ■  Ojjdt  (25) 

where  Oja  are  constant  rates  and  time  t  is  measured  from  the  beginning  of  the 
interval.  Equation  (21)  can  now  be  reduced  to  an  ordinary  differential  equation 
over  a  small  time  interval.  This  is  done  by  substituting: 


/■jr- . 


/ 


3K 

3e: 


••p  *  no  r  ?/x  3f  5?  *p  _  3f  • 

C  ^  3/2  C  2/3  fft?  ^  C  -  k  CT .  . 

•p  30^  30^  3a^7  ij 


(27) 


The  solution  to  equation  (27)  can  be  shown  to  be  of  the  following  form: 


*  At  +  k,  +  k-e 


-Xt 


(28) 


where  by  definition: 


j 


1 


j 


i 


A  = 


j  (: 


3/2  C  ^  2/3  ^ 


3f  ) 

3a. .  3a. .  / 

ij  ij  • 


X  -  3/2  C 


/l^!r 


(29) 


and  k}  and  k2  are  unknown  constants.  Time  t  is  measured  from  the  start  of 
the  interval.  The  constants  ki  and  kj  are  evaluated  from  the  conditions  at 
the  start  of  the  interval: 


at  t  =  0 


ep  «  eP 

o 


•p  »p 

«  E- 


(30) 


Using  equations  (30)  to  evaluate  k^  and  k2  gives: 

kl  *  £o  -  X  <A  - 

k2  ■  X  (A  •  £0>  <3‘) 

Substituting  into  equation  (28)  for  kj  and  k2  results: 

eP  *  At  ♦  eJ  ’  J  (A  -  eJ)U  "  e‘Xt)  (32) 
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Differentiating  equation  (32)  gives  the  rate  of  change: 

ep  *  A  -  (A  -  ep)e"Xt  (33) 

Consider  now  a  time  interval  t«0  to  t»At  and  use  equation  (33)  to  obtain  the 
change  of  the  effective  plastic  strain  rate: 


de 


P 


-  (A  -  lp)  (1  -  e_XAt)  (34) 

Returning  to  the  incremental  stress -strain  relation  and  substituting  first 

for  the  rate  term  from  equation  (19)  and  then  expressing  dep  from  equation  (34) 
gives  the  following: 


do. 


Aijkldekl 


♦  DE 


3f  3K 


ijrs  3ors  3-p 


(A 


ep)(l  -  e*^) 


(35) 


Equation  (35)  is  now  a  suitable  incremental  relation  over  a  load  time  step 
At  which  gives  the  change  in  stress  in  terms  of  change  in  total  strain  and  in 
terms  of  parameters  which  can  be  calculated  from  previous  time  step.  The 
second  term  on  the  right  hand  side  of  equation  (35)  represents  the  rate  effects 
and,  in  the  finite  element  model,  it  will  contribute  to  the  body  force. 
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III.  NUMERICAL  CALCULATIONS 


The  analytical  development  of  the  previous  section  has  been  incor¬ 
porated  into  the  elastic-plastic  version  of  the  SANX  computer  code.3 
The  basic  arrangement  of  the  SANX  code  has  been  retained.  The  main  changes 
to  the  program  involve  changes  in  the  Subroutine  ELPLSS  which  assembles 
the  incremental  plastic  stress-strain  relations.  In  the  new  version  the 
incremental  stress-strain  relations  are  based  on  the  equation  (35).  The 
basic  input  procedure  for  the  new  SANX  is  the  same  as  the  original  code 
except  for  the  first  input  card.  This  card  has  been  modified  to  input 
the  rate  dependent  material  parameters  b,  e  defined  in  equation  (2),  and 
the  time  interval  At  needed  for  the  time  dependent  incremental  solution. 

The  first  input  card  in  the  original  elastic-plastic  SANX  was: 

Card  1  (Original) 

Format  (2110) 

Columns  1-10  NTOTS 

Number  of  segments  (8  maximum) 

11-20  NOLINC 

Number  of  load  increments 


The  new  input  card  is: 

Card  1  (New) 

Format  (2U0.3E12.6) 

Columns  1-10  NTOTS 

11-20  NOLINC 

21-32  DELTIM 

Time  increment  At 

33-44  BVR 

Material  parameter  b 

45-56  EVR 

Material  parameter 

In  order  to  check  the  new  program  constant  stress  rate  was  applied 
to  sample  examples  which  simulates  uniaxial  loading  and  the  results  were 
compared  to  those  obtained  from  uniaxial  solution  of  Reference  1.  The 
uniaxial  analysis  of  Reference  1  is  in  two  parts.  The  first  part  is  an 
incremental  uniaxial  solution  and  the  second  is  an  exact  solution  for 
constant  stress  or  strain  rate  loading.  The  incremental  solution  of 
Reference  1  was  first  compared  to  the  exact  solution  and  it  was  found 
that  modifications  were  necessary  to  the  incremental  solution  to  make  it 
agree  with  the  exact  formulation.  The  comparison  of  the  results  from  the 
three-dimensional  SANX  program  are  made  to  the  modified  incremental, 
uniaxial  formulation. 
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The  first  comparison  was  made  between  the  uniaxial  incremental 
solution  and  the  uniaxial  exact  solution  from  Reference  1.  The  purpose 
of  this  comparison  was  to  establish  the  effect  on  the  accuracy  of  the 
size  of  the  load  steps  at  various  stress  rates.  Figure  1  shows  these 
results  where  the  incremental  and  the  exact  solutions  are  shown  for  a 
wide  range  of  stress  rates.  The  results  are  presented  both  in  Metric 
and  British  units.  The  Metric  units  are  given  in  parentheses  directly 
under  the  corresponding  British  units.  Figure  1  presents  results  for 
two  load  steps  of  2x10*  psi  (13.7  GPa)  and  8xl03  psi  (55.1  GPa).  The 
time  step  At  was  adjusted  to  give  the  desired  stress  rate.  These  results 
are  for  the  following  material  parameters: 

rate-dependent  yield  parameters; 

b  *  3.67  x  10"^ 

e  =3.0  x  10'^ 

0 

elastic  modulus; 

E  =  16.8  x  106  psi  (115.7  GPa) 
strain  hardening  parameter; 

C  =  0.259  x  106 
yield  stress; 

cy  =  133  x  103  psi  (0.916  GPa) 

Comparing  solutions  in  Figure  1  it  can  be  seen  that  relatively  good  agree¬ 
ment  exists  between  exact  and  incremental  solutions.  This  is  especially 
true  for  the  smaller  load  step  of  2x10^  psi  (13.7  GPa).  It  is  expected 
that  smaller  load  steps  will  give  more  accurate  results. 

The  next  comparison  involved  using  the  new  three-dimensional  computer 
program  to  analyze  a  uniaxial  situation  composed  of  a  cylindrical  body 
subject  to  axial  load.  The  results  are  compared  to  the  exact  uniaxial 
solution  of  Reference  1.  These  results  are  presented  in  Figure  2  for  the 
two  different  load  step  sizes  used  in  Figure  1.  As  in  Figure  1,  the  resuLts 
are  presented  in  two  different  systems  of  units.  It  is  expected  that  the 
results  in  Figure  2  should  duplicate  results  in  Figure  1  if  the  three- 
dimensional  code  is  working  properly.  It  can  be  seen  that  the  results 
between  uniaxial  incremental  solution  and  the  finite  element  program  are 
almost  identical. 
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APPENDIX  A 


Computer  Listing  for  the  Program 
SANXVR  for  the  Analysis  of  Nonaxi- 
syametric  Configuration  with  Rate 
Dependent  Yield  Criterion 


c*  * 
c 

c*  * 


PROGRAM  SANXVRC INPUT * OUTPUT * TAPES- INPUT , TAPES -QUTP 
1  TAPE2 , TAPE 3  k  TAPE  15  > 

2T  APE21 • T  APF25 * T  APE26 > 

a:******#***:*****#********* 
BRLESC  FINITE  ELEMENT  STRESS  ANALYSIS  OF  AXI3YMMETI 
PLANE  STRAIN*  AND  PLANE  STRESS  SOLIDS  UITH  CRTH07RC 


•  OUTPUT  iTAPEl  * 


TEMPERATURE-DEPENDENT  MATERIAL  PROPERTIES 
#$#**#***#*********#**#%  *  *  i  * 

INTEGER  CODE 

COMMON/VISC/EPSDN!  12*  10,8  >»SIGVP!  6  )*DEPSR(  6»  10 * 8  >»BELTXK 
COMMON/RATE/DKPR , SIGPR  *  BVR *  EUR , PSRATE!  10, 8 ) * NRATE 
COMMON/INCR/NOL INC  *  NOL  *  INERT , NUMMAT  *  SIGTOtt  12  k  4  >  G  ) 

1  *EF'STOT!  12 »  4k8  ) 

COMMON/PLAS/ALFA! 6*  4*8  ),SIGYLD! 7*6*3  )k IFGPL(  4k3) 

COMMON/BOCON/NRDF*NREQ! 18  >,URES! 18 ) 

COMMON/NPDATA/R(  10  )»CODE<  10  ),XR!10  )»Z!  10  >*XZ!  10  ) 

1NPNUM!  4 »  8 ) » T!  10  ),XT!  10  ) 

COMMON/ ARG/RRR<  5  )  k  ZZZ(  5  )  *  RR<  4  )  k  ZZC  4  )  k  S!  15  k  15  )  k  FT  15  ) .  TT!  6 )  » 
1H(  6*15  )»CRZ!  6*6)*XI!  10  )»ANGLE!  4  ),SIG!  18)kEPS(  18)kN 
COMMON/ELDATA/BETA!  10  ),EPR!10  )»PR!4  ),SH!4  )»IX!8 


) » JT( 4 


*  5 )» 
)* 


VCL ,  NUMNP *  NUMEL  *  NUHPC » NUM3C » 


IIP!  4  )*  UP!  4  )*  IS!  4  )» JS!  4  )*ALF‘HA(  10  ),IT!4  )»JT<4  )» 

2ST(  4  ) 

COMMON/BASIC/ACELZ* ANGVEL*ANGACC, TREE  ,  VCL»NUMNP*NUMEL*NUMPC  »NUM3C» 
1NUMST 

COMMON/NXMESH/THETAN!  8  )* NPC!  8*8 ) 

C0MM0N/ANS1/NUMELS(  8 )*NUMNPS!  8 ) 

COMMON/NXDATA/NTP * NTS  »NT0TS  *  GTS1G(  24,24,8) 

COMMON/NONAXI/Sl!  30,30  ),P1<  30  ), THETA , BS1(  6,30  ) 

COMMON/SOL VE/X! 888  ),Y!  888  )»TEM!  888  >,NUHTC* MBAND 

COMMON/TD/IMIN<  100  ) ,  IMAX!  100  ) ,  JMIN!  25  ) ,  JMAX(  25  ) , MAXI , MAX JVNMTL*  NBC 

COMMON/CONVRG/IDONE 

COMMON/PL ANE/NPP 

COMMON/RESULT/BS!  6 , 15  ) , D<  6 , 6  > , C! 6  *  6  ) , AR , BBC  6 , 9 ) ,  CNS(  6,6) 
COMMON/MATP/RO(  6  ),E<  12, 16*6  ),EE>’  16  ),AOFTSC  6  ) 

COMMON/NXSOLV/SKG!  36  »24>*FTG. l32),FT0T( 132),IT0T 
DIMENSION  TITLE( 20  ) 

*#####**###  *####**  it!#:***  *********** 


READ  AND  WRITE  CONTROL  INFORMATION 
He********************; 
RE AD( 5,3000)  NTOTS , NOL I NC , DELTIM , BUR ,  I 
WRITE! 6,3017  )  BVR, EUR 
017  FORMAT! 1H  »"  BVR  =  " *E12.4»"  EVR  =  ” 
DO  150  1=1* NTOTS 
150  READ! 5*3001  )  THETAN! I > 

DO  152  1=  1, NTOTS 

152  READ!  5*3002  )  (  NPC!  I ,  J  ),  J  =  1 , 8  ) 

3000  FORMAT! 2110, 3E12, 6  ) 

3001  FORMAT! F10. 5,110  ) 

3002  FORMAT! 8110) 

REWIND  15 
REWIND  26 
REWIND  21 


.2,4  ) 


3010 


3011 


REWIND  26 
REWIND  21 
REWIND  25 
WRITE! 6,3010 ) 

FORMAT! "1" *n SEGMENT  DATA  FOR 
WRITE! 6*3011  )  NTOTS  ,NOLIMC 
FORMAT!"  NUMBER  OF  TOTAL 


OBLEM"  ) 


NONAX I SYMMETR I C  PRQBLE 
, DELTIM 

SEGMENTS  =",I5,//» 


1  NUMBER  OF  LOAD  INCREMENTS  =" ,15,//, 
TIME  INCREMENT  =",E15.8) 
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DO  153  1=1 .NTCTS 
WRITEt  6,3012  )  I»THETAN<  I) 

3012  FORMATC1  “  ,///,"  SEGMENT  TYPE  =",I5/»B  THETA  =  "  ,F10.5) 

153  CONTINUE 

DO  154  I=1»NT0TS 

154  WRITE<  6»3014)I»(  NPC<  I » J  )» J=1 ,8  ) 

3014  FORMATC1  "»" CONNECTING  NODES  FOR  SEGMENT11  *  I5»  "  ARE" ,815) 

DO  910  NOL=lf NOLINC 
WRITE(  6,2030  )  NOL 
REWIND  15 

DO  950  NTP  =  1»NT0TS 
THETA=  THETAN(NTP)  757. 295730 

IF< NOL  »NE ♦ 1 >G0  TO  525 

50  RE AIK  5 , 1000  >TITLE , NNL A , NUMTC , NUMMAT  »  NUMPC » NUMSC , NUMST » TREF 
1 , 1 NERT , NL INC » INC I , I NCF , I PLOT 

URI TE(  6 , 2000  )T ITLE , NNLA  ,  NUMTC ,  NUMMAT » NUMPC , NUMSC , NUMST  # TREF , INERT , 
1NLINC 

URITEt 15 )NUMTC, NUMMAT, NUMPC, NUMSC, TREF, INERT,  INCI » INCF 
NPP=0 

C*  *********************************** 
C  GENERATE  FINITE  ELEMENT  MESH 

C*  ***************************  ******** 
100  CALL  MESH 

DO  155  1=1 . NUMEL 
IFGPL< I, NTP >=0 
PSRATE< I, NTP  )=0 .0 
DO  155  J=1 .12 
SIGTOT(  J.I.NTP  )=0.0 
EPSTOT< J.I.NTP  )=0.0 
ALFA< J.I.NTP  )=0.0 
EPSDN< J.I.NTP >=0.0 

155  CONTINUE 

WRITE<  15  X  R<  I  ),  1=1  .NUMNP  ) 

URITE<  15  X  Z<  I  ), 1  =  1, NUMNP) 

NUMELS(NTP)  =  NUMEL 
NUMNPS(NTP)  =  NUMNP 
IF  (IPLOT.EG.l)  CALL  MPLOT 

C*  *********************************** 
C  READ  AND  WRITE  TEMPERATURE  DATA 

C*  ************************  *********** 
103  IF< NUMTC. EB.O)  GO  TO  440 

IF< NUMTC. GT.O)  READ(5,1001)  <  X<  I )» Y<  I  ),TEM(  I )» 1=1  ,NUMTC ) 

IF< NUMTC. EO. -2)  CALL  TEM2< NUMNP ) 

IF<  NUMTC.EG.-2  )  GO  TO  440 

MPRINT=0 

DO  210  1=1, NUMTC 

IF< MPRINT.NE.O)  GO  TO  200 

WRITE<  6,2001  ) 

MPRINT=59 

200  MPRINT=MPRINT-1 

210  WRITE<  6,2002)  X<  I  >, Y<  I  ),TEN<  I  > 

MPRINT=0 

DO  230  N=l, NUMNP 

IF(  MPRINT.NE.O)  GO  TO  220 

WRITE<  6,2003) 

MPRINT=59 

220  MPRINT=MPRINT-1 

CALL  TEMP(  R<  N  ),Z<  N  ),T<  N  ) ) 

230  WRITE(  6,2004)  N,R<  N  ),Z<  N  )»T<  N  ) 

440  MPRINT=0 
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DO  460  N=1 .NUMEL 
IFC  MPRINT .NE.O  >  GO  TO  450 
„  WRITEC  6.2008  ) 

MPRINT =59 

450  MPRINT=MPRINT-1 
r-  II  =  IXCN.l) 

J.J=IXCN.2) 

KK=IXC  N » 3  ) 

LL=IX<  N.4  ) 

C 

C  TEM  IS  TEMPORARY  STORAGE  FOR  ELEMENT  TEMPERATURES 

C 

TEM<  N  )=(  T<  II  )+TC  JJ  HTC  KK  )+T(  LL  )  >74.00 
460  WRITEC  6*2009  )  N.C  IX<  N»I  ) » 1=1 .5  J.BETAC  N  J.ALPHAC  N  ).TEM(  N  ) 

WRITE<  15  )( <  IX<  I»J)»J=1»5)»  1  =  1 » NUMEL  ) 

WRITEC  15  )<  BETA  C I ), 1=1 .NUMEL > 

WRITEC  15  )<  ALPHA(  I )» 1=1 .NUMEL > 

WRITEC  15  )(  TEM  ( I  )» 1=1 .NUMEL  ) 

DO  470  K=l. NUMEL 
470  T<  K  )=TEM<  K  ) 

C*  ****************  ******************* 

C  READ  AND  WRITE  MATERIAL  PROPERTIES 

C*  *********************************** 
500  CONTINUE 

DO  510  M=1 .NUMMAT 

REAIK  5.1004)  MTYPE .  C  NT  > ROC  MTYPE  ) » AOFTSC  MTYPE  ) ) 

WRITEC  6 » 2010)  MTYPE . NT . ROC MTYPE  ) 

REAIK  5.  1005  )C  C  EC  I » J.MTYPE  )» J=1 » 14 >» 1=1 »NT  ) 

REAIK  5.1005  )C  SIGYLDC  I .MTYPE.NTP  >.1=1.7) 

IFC  AOFTSC  MTYPE  )«NE»1 *  )  WRITEC  6.2011  )C  C  EC  I .J.MTYPE  >, J=l,13>.  1=1. NT  ) 
IFC  AOFTSC  MTYPE  )»EQ.l »  )  WRITEC  6. 2012  >C  C  EC  I .  J, MTYPE  > » J=1 . 13  > » 1=1 » NT  > 
WRITEC  6.3015  )C  SIGYLDC  I .MTYPE.NTP  )» 1=1 .7  ) 

3015  FORMATC 1H  ."YIELD  STRESSES  ARE  t"»7 
11H  » " Y 1 1  =  " »E15 .7/ 

21H  ."Y22  =  B » E15.7/ 

31H  »"  Y33  =  " »E15.7/ 

41H  » "  Y12  =  " .E15.7/ 

51H  » "  Y13  =  " »E15 .7/ 

61H  » "  Y13  =  " .E15.7/ 

TIM  ■  »  P  —  7* 

WRITEC 15  )MTYPE. NT »  ROC MTYPE  ) 

WRITEC  15  )C  C  EC  I .J.MTYPE  >»J=1»14  )»I=1»NT> 

DO  510  I=NT .12 
DO  510  J=1 » 16 

510  EC  I » J  »  MTYPE  >=EC  NT » J » MTYPE  ) 

GO  TO  526 

525  CALL  DATA 

526  CONTINUE 

DO  900  NL=1»NLINC 
ACELZ=0.00 
ANGVEL=0.00 
ANGACC=0.00 

IFC INERT  .EQ.  0 )  GO  TO  511 

IFC  NL  .NE .  1  » AND ♦  INCI  .EQ.  0 )  GO  TO  511 

C  READ  AND  WRITE  DYNAMIC  FORCES 

*********************** 

RE ADC  5.1030  )  ACELZ.  ANGVEL.  ANGACC 
WRITEC 6.2031)  ACELZ.  ANGVEL.  ANGACC 

511  CONTINUE 


*  * 


RITE  PRESSURE  AND  SHEAR  BOUNDARY  CONDITIONS 


#  If.  ^  n'  »V  *  n'  *V  V- 

C  READ  AND  U 

C*  **************  *  * 

IF<  NL  . NE ♦  1  .AND.  INCF  .EQ.  0)  GO  TO  700 
600  IF<  NUHPC.EO.G )  GO  TO  630 
•’  MPRINT=0 

DO  620  L=1  * NUMPC 

IF( HPRINT ♦ ME . 0  )  GO  TO  610 

WRITE*  6*2013 ) 


*  *  * 
r. » ; t  ■ 

L)1L>. 

*  *  *  'i'-  *  *  n-  *  *  *  *  * 


*  *  *  * 


MPRINT-53 

610  MPR I NT =MPRINT- 1 

READ*  5*1006  )  IP*  L  >»  JP<  L  ) *  PR(  L  ) 
620  WRITE*  6*2014  )  IP*  L  ) *  JP*  L  ) > PR*  L  ) 
630  IF< NUMSC.EG.O  )  GO  TO  701 
MPRINT=0 


DO  650  L=1 »NUMSC 

IF( MPRINT.NE.O)  GO  TO  640 

WRITE* 6*2015) 

MPRINT=58 

640  MPRINT=MPRINT-1 

READ*  5*1006  )  IS*  L  )» JS*  L  ),SH*  L  ) 

650  WRITE*  6*2014)  IS<  L  )*JS<  L  ),SH<  L  ) 

701  IF*  NUMST .EQ.O  )  GO  TO  700 
MPRINT=0 

DO  680  L=l* NUMST 

IF* MPRINT.NE.O)  GO  TO  670 

WRITE*  6*2025) 

MPRINT=58 

670  MPRINT=MPRINT-1 

READ*  5, 1006)  IT*  L  )*  JT<  L  )»ST(  L  ) 

680  WRITE*  6*2014  )IT*  L  )*  JT*  L  )*ST*  L  ) 

c*  *********************************** 

C  DETERMINE  BANDWIDTH*  INITIALIZE  ELASTIC-PLASTIC  RATIO, 

C  AND  CONVERT  BETA  FROM  DEGREES  TO  RADIANS 

C*  *********************************** 
700  J=0 

DO  710  N=1 *NUMEL 
IX<  N,5  )=IABS<  IX<  N * 5  ) ) 

DO  710  1=1,4 
HO  710  L=l»4 

KK=IABS(  I X<  N  *  I  )- 1 X<  N  *  L  ) ) 

IF( KK.GE.J  )  J=KK 
710  CONTINUE 

MBAND=3*J+3 
IF( NL.GT.l  )  GO  TO  721 
DO  720  N=1 *NUMEL 
EF'R<  N  )=1 . 

ALPHA<  N  )=ALPHA<  N  )/57 . 295780 

720  BETA(N)=BETA(N)/57. 295780 

721  CONTINUE 

C*  **********************  ************* 
C  SOLVE  NONLINEAR  PROBLEM  BY  SUCCESSIVE  APPROXIMATIONS 

C*  ***************************  ******** 
no  800  NNN=1 »NNLA 
C 

C  FORM  STIFFNESS  MATRIX 

C 


CALL  STIFF 
C 

C  SOLVE  FOR  niSPLACEMENTS 
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c 


CALL  SOLV 


C 

C  COMPUTE  STRESSES 

C 

*  CALL  STRESS 

CALL  STORE 

IF< I DONE ♦  NE  •  1  )  GO  TO  300 

799  NITER=NNN 

IF(  IDONE.EQ. 1  )  GO  TO  310 

800  CONTINUE 

810  IF* IDONE.EQ. 1  )  WRITE* 6*2016  )  NITER 
IF< IDONE.NE. 1  )  WRITE* 6*2017)  NITER 
900  CONTINUE 
950  CONTINUE 


IT0T=24+12#< NTOTS-1  )  NEWDYN 

IF*NOL.NE.l)  GO  TO  83  NEWDYN 

C******%********%**%*%%%***%*%*****M*M*%K*MM%t***t***%****%MM*t*  NEWDYN 
C  INITIALIZE  PREVIOUS  HISTORY  TOTAL  DISPLACEMENTS  NEWDYN 

DO  89  1=1  * ITOT  NEWDYN 

FTOT*I)=0.00  NEWDYN 

89  CONTINUE  NEWDYN 

88  CONTINUE  NEWDYN 


CALL  ASEMBL 
CALL  ANSWER 
910  CONTINUE 

1000  FORMAT*  20A4/6I5»F5. 0*515  ) 

1001  FORMAT* 3F10.0  ) 

1004  FORMAT  *2I5»2F10.0) 

1005  FORMAT* 7F10.0  ) 

1006  FORMAT  <2I5»F10.0) 

1030  FORMAT* 3F 10.0  ) 

2000  FORMAT  * 2H1  *20A4/ 

1  33H0  NUMBER  OF  APPROXIMATIONS - 14/ 

2  33H0  NUMBER  OF  TEMPERATURE  CARDS  — 14/ 

3  33H0  NUMBER  OF  MATERIALS - 14/ 

4  33H0  NUMBER  OF  PRESSURE  CARDS - 14/ 

5  33H0  NUMBER  OF  SHEAR  CARDS - 14/ 

6  33H0  NUMBER  OF  TORSION  CARDS - 14/ 

7  33H0  REFERENCE  TEMPERATURE - El 2. 4/ 

8  33H0  NUMBER  OF  INERTIA  CARDS - 14/ 

9  33H0  NUMBER  OF  LOAD  INCREMENTS - 14/) 

2001  FORMAT  * 1H1 »13X* 1HR* 14X * 1HZ* 14X  * 1HT  ) 

2002  FORMAT  (3F15.3) 

2003  FORMAT  < 35H1  NR  Z  T) 

2004  FORMAT  ( I5*2F10.4*F10.0  ) 

2008  FORMAT  <  74H1  EL  I  J  K  L  MATERIAL  ANGLE  BETA  ANGLE  A 
1LPHA  TEMPERATURE  ) 

2009  FORMAT  *  15* 414 *  18 *F1 1 . 1 , 2F13 . 3  ) 

2010  FORMAT  < 1H1 » "MATERIAL  IDENTIFICATION  NUMBER  =",12/ 


1 1 H 
21H 


’■NO*  OF  MATERIAL  TEMPERATURE  CARDS  ="»I2/ 
1  “MASS  DENSITY  =  "»E15.7) 

*  2011  FORMAT  <  1H  » " TEMPERATURE  ="»E15.7/ 

11H  » "MODULUS  OF  ELASTICITY-EN  =",E15.77 
OF  ELASTICITY-ES  ="»E15.7/ 

OF  ELASTICITY-ET  =” ,E15.7/ 
RATIO-NUNS  = " » E 1 5 . 7 / 

RATIO-NUNT  =",E15.7/ 


21H 

31H 

41H 

51H 

61H 

71H 

81H 

91H 

11H 

21H 

31H 


2012  FORMAT  (  1H  » 


11H 

21H 

31H 

41H 

51H 

61H 

71H 

81H 

91H 

11H 

21H 

31H 


"MODULUS 
"MODULUS 
"POISSON 
"POISSON 
"POISSON  RATIQ-NUST  =" fE15.7/ 

"SHEAR  MODULUS-GNS  =”,E15.7/ 

"SHEAR  MODULUS-GST  =M»E15.77 
"SHEAR  MODULUS-GTN  =" »E15.7/ 
"COEFFICIENT  OF  THERMAL  EXPANSION-AN 
"COEFFICIENT  OF  THERMAL  EXPANSION-AS 
"COEFFICIENT  OF  THERHAL  EXPANSION-AT 


=" »E15.7/ 
=" »E15.7/ 
=" f E15*7/  ) 


TEMPERATURE  =" .E15.7/ 
"MODULUS  OF  ELASTICITY-EN  =".E15.7/ 
"MODULUS  OF  ELASTICITY-ES  =",E15.7/ 
"MODULUS  OF  ELASTICITY-ET  =%E15.7/ 
"POISSON  RATIO-NUNS  =" >E15.7/ 
"POISSON  RATIO-NUNT  ="»E15.7/ 
"POISSON  RATIO-NUST  ="»E15.7/ 

"SHEAR  MODULUS-GNS  ="»E15.77 
"SHEAR  MODULUS-GST  =" .E15.7/ 

"SHEAR  MODULUS-GTN  ="»E15*7/ 

"FREE  THERMAL  STRAIN-FN  =",E15.7/ 
"FREE  THERMAL  STRAIN-FS  =",E15.7/ 
"FREE  THERMAL  STRAIN-FT  ="»E15.7/> 


2013 

2014 

2015 

2016 
2017 

2024 

2025 

2030 

2031 


<  30H  PRESSURE 
( 2I5»F10. 1 ) 


BOUNDARY  C0NDITIQNS/20H 


J  PRESSURE ) 


(  27H 

<  26H 

<  33H 
(43H0 


SHEAR  BOUNDARY  CGNBITI0NS/17H 
THE  SYSTEM  CONVERGED  IN  I2»11H 
THE  SYSTEM  DID  NOT  CONVERGE  IN 


I  J  SHEAR ) 
ITERATIONS) 

ITERATIONS) 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 

FORMAT  (43H0  THE  AXISYMMETRIC  OPTION  NAS  BEEN  SELECTED) 

FORMAT<  30H  TORSION  BOUNDARY  CONBITIONS/17H  I  J  SHEAR) 
FORMAT*  1H  »45Xf ” ************  LOAD  STEP  *************  ="iI4) 
FORMAT(  1H0  , “AXIAL  ACCELERATION  =" rE12.4/ 

11H0  , "ANGULAR  VELOCITY  ="»E12.4/ 


21H0  , "ANGULAR  ACCELERATION^’  f E12.4 ) 

920  STOP 
END 

SUBROUTINE  ANGLE  <  R*Z>RC>ZC»ANG ) 

C  FIND  ANGLE  OF  INCLINATION  BETWEEN  0  AND 
C*  ********************** 
PI-3 . 1415927 
I«1=(Z-ZC) 

D2=( R-RC  ) 

IF(  ABS(  R-RC  )»GT .  1  *E-8  )  GO  TO  100 
ANG-PI/2. 

IF< Dl.GT.l.E-B)  RETURN 
ANG=-ANG 


2*PI 
*  * 


*********** 


RETURN 

C*  *********************************** 
C  ALLOW  CIRCLE  TO  CROSS  AXIS 

C*  **************************  ********* 
100  ANG=ATAN2<  Dl» D2  ) 

RETURN 

END 

SUBROUTINE  ANSWER 
INTEGER  CODE 


24 


COMMON/VISC/EPSDN<  12*  10,3  >»SIGVP<  6  )*DEPSR(  6,10,3>»r!ELTIM 
COMMON/RESULT/BS<  6 , 15  ) * B<  6  *  6  >  »  C(  6 * 6  )  *  AR » BBC  6 , 9  ) » CNSC  6,6 ) 
C0MM0N/PLAS/ALFA<6,  4# 8  )*SIGYLEK 7*6, 8>, IFGPL(  4*8  > 
COMMON/INCR/NOLINC » NOL » INERT » NUMMAT *SIGTQT<  12.  4*0) 

1  ,  EPSTOT< 12,  4,8  ) 

COMMON/ELDATA/BETA(  10  >,EPR(10  >,PR<4  >,SH<4  >,IXC8  ,5)» 

1  IP<  4 '  ),JP<  4  )*  IS<  4  >*JS<4  )*  ALPMAC 10  >,IT(4  ),JT<4  >, 

2ST<  4  ) 

COMMON/ARG/RRR<  5 ),ZZZ<  5  >,RR<  4  >»ZZ<  4  )»S(  15,  15  >*  P(  15  )»TT(  6  ), 

1H(  6*15  ),CRZ<  6*6 ),XI(  10  >*  ANGLE<  4  )»SIG<  18)*EPS<  18>*N 
COMMON/NXSOLV/SKG( 36  , 24  >  *  FTG< 132  > , FTOT(  132  ) » ITOT 
COMMON/ ANS2/  UT1<24)»  G<24*24>,  GR1(  24*24  ) »  BUMM(  24*24  > 
C0MM0N/ANS1/NUMELS<  8  )»NUMNPS< 8  ) 

COMMON/NONAXI/SK  30,30  >*P1<  30  >»THETA*BS1<  6  ,30  ) 

COMMON/NXDATA/NTP , NTS , NTOTS * GTS1G(  24 *  24,8  > 

COMMON/NXMESH/THETAN< 8  )  » NPC<  8.8) 

COMMON/ARGl/SIGK  18>,EPS1<  18>,BEPSP<  12  >*CEPSP<  6,6 ) 

COMMON/SOLVE/B<  72>»A<  72*36 >*NUMBLK*MBANB 
DIMENSION  UT<  24  >»UC1<  24  )*UC<  24  >*R1(  24*24  ) 

REWIND  25 
REWIND  26 
REWIND  21 
K0LD=1 

DO  100  K=l, NTOTS 

NS=K 

KNEW=K 

NUMNP  =  NUMNPS<  K  > 

NUMNP3  =  3*NUMNP 
NUMEL  =  NUMELS(K) 

K20  =  21 

READ(  26  )  <  B(  I  )» 1=1 » NUMNP3  ) 

REAIK  26  >  ( (  IX<  I  *0  > ,  J=1 ,5  >»  1=1 » NUMEL  ) 

WRITEt 6,1200)  K 

REAIK  25  )( <  Rl<  I ,  J  ) ,  J= 1 » 24  ) » 1= 1  * 24  ) 

DO  110  KK=1 *4 
NP1  =  NPC( NS»KK  ) 

NP2  =  NPC<NS,KK+4> 

DO  110  1=1,3 

UC<  3*<  KK-1  >+I  )  =  B<3*NPl-3+I> 

UC(  3*(  KK-1  )+I+12  )  =  B<3*NP2-3+I) 

110  CONTINUE 

DO  115  KK=1 ,24 

115  UT<KK)  =  FTG(  KNF(  NS-1  )#12  ) 

WRITE(  6,900) 

900  FORMAT< *  EL  SIGMAR  SIGMAZ  SIGMAC  SIGMARZ  SIGMAZC” 

1  ,"  SIGMACR  SIGMAN  SIGMAS  SIGMAT  SIGMANS" , 

2  "  SIGMATN" ,/"  EPSR  EPSZ  EPSC  EPSRZ  ", 

3  "EPSZC  EPSCR  EPSN  EPSS  EPST  EPSNS  " , 

4  "  EPSST  EPSTN"  ) 

IF<  KOLB. EG. KNEW  )  REWIND  21 
IF(KOLD.NE.KNEW)  KOLD=KNEW 
DO  120  N=l, NUMEL 
MTYPE=IABS<  IX(  N,5>) 

READ<  K20  )( <  CRZ<  I ,J >, J=1 » 6 > *  1=1 ,6  ) 

READ(  K20  )<  <  BS1<  I, J),J=  1,30),  1=1,6) 

RE AD<  K20  X  <  G<  I ,  J  ) ,  J=  1 , 24  ) » 1  =  1 , 24  ) 

READ(  K20  )( <  CEPSP<  I, J  )» J=l*6  >,1=1,6 > 

READ<  K20  >(  <  CNS<  I,J),J=1,6)*I=1*6> 

READ(  K20  >(<D<  I,J),J=1,6),I  =  1*6) 

READ(  K20  )(<C<I,J),J=1,6),I=1,6) 
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DO  125  1=1,24 
DO  125  J=1 ,24 

GR1C I,  J  )  =  0.00 
DO  125  M=1 » 24 

125  GR1C  I » J  )  =  GR1C  I » J  )  !  G(  I  >  M  )4R1C  M»  J  ) 

DO  126  1=1,24 
UC1( I  )  =0.00 
UT1( I  )  =0.00 
DO  126  J=1 ,24 

UC1(  I)  =  UCKI)  4  GR1C  I,J  )*UCC  J) 

126  UTKI)  =  UT1(  I  )  4  GR1C  I ,  J  )*UT<  J  ) 

DO  130  1=1,4 
11=3*1 

JJ=3*IXC N , I  ) 

PIC  II-2)  =  B<  JJ-2) 

Pl(  II-l  )  =  BC  JJ-1  ) 

Pl<  II  )  =  BC  JJ  ) 

Pl<  1 1410)  =  B<  JJ-2  ) 

Pl< 1 141 1  )  =  B< JJ-1  ) 

PIC  1 1412  )  =  B<  JJ>  / 

130  CONTINUE 

DO  135  1=1,24 

135  Pl(  I  )  =  Pl<  I )  -UC1C  I  HUT1C  I  ) 

DO  136  1=1,3 

Pl(  1424  )=  (  PIC  I  )4P1<  143  )4P1C  146  )4P1C  149  )  >/4.00 

136  PIC  1427  )  =  C  PIC  1412  )4P1C  1415  )4P1C  1418  )4P1C  1421  ))/4. 00 

DO  140  1=1,6 
EPS1C I  )  =  0.00 
DO  140  J=1 ,30 

140  EPS1CI)  =  EPS1C  I  >4BS1C  I ,  J  )*P1C  J  ) 

DO  143  1=1,6 
EPS 1C  146  )=0.0 
DO  143  J=1 ,6 
DO  143  L=1 ,6 

143  EPS1C  146  )=EPS1C  I46)4DC  I ,  J  )*CC  J , L  )*EPS1C  L  ) 

DO  150  1=1,6 
SIG1C I  )  =  EPSBNC I,N,NS) 

SIG1C 146  )=EPSDNC 146, N, NS  ) 

SIGVPC  I  )=0.0 
DO  150  J = 1 , 6 

SIG1C  I  )=  SIG1CI)  4  CRZCI,J)*EPS1CJ) 

150  SIG1C 146  )=SIG1C 146  )4CNSC I,J  )*EPS1( J46) 

DO  151  1=1,6 

SIGVPC  I  )=SIG1C  146) 

151  CONTINUE 

DO  141  J=1 , 12 

141  EPS1CJ)  =  EP51C J  )*100.0 
DO  230  1=1,6 

DEPSPC  I  )=0,0 
230  DEPSPC 146  )=0.0 

IF  C  IFGPLCN,NS).EQ.O)  GO  TO  241 
DO  250  1=1,6 
DO  250  J=1 ,6 

250  DEPSPC 146  )=DEPSPC 146  )4CEPSPC I,J  )*EPS1C  J46  )/100. 

DO  251  1=1,6 

251  DEPSPC 146  )=DEPSPC 146  )4DEPSRC I,N,NS) 

DC  4,1  )=0,5*DC  4,1  ) 

DC  4,3)=  0.54DC  4,3  ) 

DC  1 ,6  )=2.0*DC  1,6) 

DC2»6>=  2.0*DC  2,6 ) 
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CC  1 ,4  >=  2.0&CC  1  >4  ) 

C<  2*4  )=  2.0$CC  2*4  ) 

CC  4*1  )=  0.5*CC  4*1  ) 

CC  4»2>=  0 ♦  5$CC  4*2) 

DO  160  1-1*6 
DEPSPC  I  )=0.0 
'  DO  160  J=1 *6 
DO  160  L=1 *  6 

160  DEPSP<  I  >=BEPSPC  I  )+CC  J*  I  )*B<  J*L  )*DEPSPC  L+6  ) 

URZTE(  6*  1400  )C  DEPSP<  I  )>  1=1*12) 

1400  FORMATC "  PLASTIC  STRAINS" /2X *  12E10.4) 

C  DO  233  1=1*6 

C  EPSDNC I  *  N  *  NS )=DEPSP(  1+6  )/BELTIM 
C  233  CONTINUE 
241  CONTINUE 

DO  240  1=1*12 

SIGTOT(  I *N*NS  >=SIGTOTC I *N*NS  )+SIGl(  I > 

240  EPSTOT< I »N*NS  )=EPSTQTC I ,N*NS  )+EPSl(  I ) 

WRITEC 6*1000  )  N*C  SIGTOT< I *N» NS  ). 1  =  1 > 12  ) 

WRITEC  6*1111  ) 

1111  FORMATC "  " * " SIGVP  "  ) 

WRITEC 6*1000  )  N»<  SIGVPC I  )»I  =  1»6) 

CALL  YIELD<  N*NS*MTYPE  ) 

IF< IFGPL< N  *NS  ) . EG . 1  )  WRITEC  6 »  1300 )N  *  NS 
1300  FORMATC"  " * " ELEMENT" » 15 * " OF  SEGMENT" *  15* "HAS  YIELDED") 
WRITE< 6*1100  >  <  EPSTOTC I ,N*NS >. 1  =  1 . 12  ) 

120  CONTINUE 
100  CONTINUE 

REWIND  21 
REWIND  25 
REWIND  26 

1000  FORMATC"  " *15* 12F9.0  ) 

1100  FORMATC"  "  »5X*12F9.5> 

1200  FORMATC "1" SEGMENT  TYPE" *  15 * // * "  " * "SEGMENT  NUMBER  =  "*I5) 
RETURN 
END 

SUBROUTINE  ASEMBL 

COMMON/INCR/NOLINC  *  NOL  *  INERT  *  NUMMAT  *  SIGTOTC  12  *  4*8) 

1  *EPSTOTC  12*  4*8) 

COMMON/BOCON/NRDF , NREQC 18  )*URESC 18  ) 

COMMON/GLBSEG/FIC  24*8  )*FEC  24*3  )*UCC  24*3)»SKC  24*24*3  )  - 
COMMON/NXDATA/NTP *  NTS  *  NTOTS  *  GTS 1GC  24 , 24  *  3 ) 
COMMON/NXSOLV/SKGC 36  *24  >*FTG< 132)*FT0TC 132)*IT0T 
C0MM0N/ANS2/FCC  24  ),GC  24*24  >*GR1C  24 * 24  ) * DUMMC  24*24  ) 

ITOT=  24  +  12*C  NTOTS-1 ) 

DO  10  1  =  1  * ITOT 
FTGC I  )  =  0,00 
DO  10  J  =  1*24 
10  SKGC  I*J)  =  0.00 
DO  100  M=1»NT0TS 

C  COMBINE  FI*  FE*  AND  SK*UC  INTO  A  TOTAL  FORCE  VECTOR  FC 
C  ***#**:} (H<#******4c******4******)»!*****4:**4;44i4c*!|c!k!(i*44:4**444444i!*:*^ 
DO  55  1=1*24 
FCC I  )  =  0,00 
DO  55  J=1 *  24 

55  FCC  I  )  =  FCC  I)  +  SNC  I  *J*M  )*  UCCJ*M> 

DO  60  1=1*24 

60  FCC  I  >  =  FCC  J.  )  +FEC  I  *M  )  -FICI*M) 
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c 

c 


r 


NOW  FILL  GLOBAL  FORCE  AND  STIFFNESS  MATRICES 

-f  '■'^V'i  'l  V 't-  $  v  t  '(■  V  h  T-  $  a-  -r.  n- ''r  *v  n- 

DO  70  1=1*24 


II  =  I+<  M-l  >$12 
FTG!  II  >  =  FTG!  ID  +  FC!  I  > 

DO  70  J=I,24 

SKG(  I1»J+1-I  >  =  SKG  ( II » J+l-I )  +  SK<  I ,  J  *  M  ) 
70  CONTINUE 
100  CONTINUE 

IF  (NOL.NE.l >  GO  TO  80 


C  READ  THE  TOTAL  NUMBER  OF  RESTRAINED  DEGREES  OF  FREEDOM 
REAIK  5*  1200  )  NRDF 
WRITE<  6* 1255  )  NRDF 

C  IMPOSE  BOUNDARY  CONDITIONS  ON  RESTRAINED  B-O-F 
DO  150  NBC=1 *NRDF 

C  READ  THE  EQUATION  NUMBER  AND  THE  IMPOSED  BOUNDARY  CONDITION 
REAIK  5 * 1 250  >  NREG<  NBC  )*URES!  NBC  ) 

URITE<  6*1260  >NREQ<  NBC  )  *  URES<  NBC  ) 

150  CONTINUE 
80  CONTINUE 

DO  160  NBC=1 »NRDF 

160  CALL  XMODFY!  URES!  NBC ) *  NREQ<  NBC  > ) 

1200  FORMAT< 15 ) 

1250  FORMAT!  I5*F10.0> 

1255  FORMAT! 1H1, "NUMBER  OF  RESTRAINED  DEGREES  OF  FREEDOM  ="*110/ 
1  "  EQUATION  NUMBER  VALUE  " > 

1260  FORMAT  < "  " *5X*I5*5X*F10.2 > 

CALL  XSOLVE 
WRITE!  6*1050) 

WRITE!  6,1100  )!  FTG!  I  >*I  =  1*IT0T> 

DO  200  1=1 , ITOT 
.  200  FTOT!  I  )=FTOT!  I  )+FTG!  I  ) 

WRITE! 6,1051 ) 

WRITE!  6*  1100  >!  FTOT!  I  )*  1  =  1,  ITOT) 

1050  FORMAT!  "1"," INCREMENTAL  DISPLACEMENTS  AT  CONNECTING  NODES"/ 
1  18X,2HUR* 18X , 2HUZ, 18X* 2HUT  ) 


1 1 00  FORMAT! “  " , 3E20 . 7  ) 

1051  FORMAT! "I" TOTAL  DISPLACEMENTS  AT  CONNECTING  NODES"/ 
1  18X* 2HUR » 18X *2HUZ* 18X * 2HUT  ) 

RETURN 

END 

SUBROUTINE  CIRCLE!  ANG1*DELPHI»RSTRT*ZSTRT*RC,ZC*I, J ) 


C* 

C 

C* 


INTEGER  CODE 


COMMON/TD/IHIN(  .100  >,  IMAX!  100  ),  JMIN!  25  ) » JMAXi  25  )»MAX2 ,  HkXJ ,  NM 


COMMON/NPDATA/R!  10  )» CODE!  10  )*XR(10  >*Z!10  )*XZ(10  )» 

1NPNUM!  4,  8  )*T!  10  )»XT!  10  ) 

DIMENSION  AR<  4*  8>*AZ!  4*  8) 

EQUIVALENCE  !  R!  1 )» AR  )»(  Z(  1  )*AZ  ) 

♦  ♦at*************************;******* 
FIND  INTERSECTION  OF  LINE  AND  CIRCLE  =  NEW  R  AND  Z 
*********************************** 
ANG1=ANG1+DELPHI 

RR=SQRT< ! RSTRT-RC  )**2+!  ZSTRT-ZC  >**2 ) 

AR!  I  * J  >=RC+RR*COS<  ANG1  ) 


AZ!  I*J  )=ZC+RR*SIN<  ANG1  ) 

RETURN 

END 

SUBROUTINE  DATA 
INTEGER  CODE 

COMMON/INCR/NOLINC»NOL» INERT  *NUMMAT*SIGTOT<  12*  4,0) 
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1  *EPST0T<12»  4,8) 

CO.MMON/NPBATA/R<  10  >, CODE*  10  ),XR<10  ),Z(  10  >,XZ(10  >, 

1  NPNUMt  4,  8  ),T<  10  ),XT<10  ) 

COMMON/BASI C/ ACELZ ,  ANGVEL ,  ANGACC ,  TREF ,  VOL ,  NUMNP »  NUMEL  ,  NUMPC  *  NUMS 
1NUMST 

COMMON/ELDATA/BETA(  10  >,EPR<  10  )>PR<  4  ),S!«4  >,IX<3  ,5)> 

1IP<4  )» JP<  4  )» IS<  4  )*  JS(  4  )»ALPHA<  10  >*IT<4  ),JT(4  >, 

2ST(  4  ) 

COMMON/NXBATA/NTP , NTS , NTOTS » GTS1G<  24 » 24 ,8 ) 

COMMON/MATP/RO<  6>,E(  12,16,6)>EE(  16),A0FTS(6> 

COMHON/SQLV£/X<  888  ) , Y<  388  ) »  TEM<  888  ) ,  NUMTC,M BAK'D 
COMMON/TD/IMIN<  100  >,IMAX<  100  >,  JMIN<  25  ),  JMAX(  25), MAXI, 

1  HAXJ,NMTL,NBC 

RE  AIK  15  >NUMTC ,  NUMMAT ,  NUMPC ,  NUMSC ,  TREF ,  INERT , 

1  INCIfINCF 

REAB<  15  >NBC,NMTL 

REAIK  15  >NUMEL, NUMNP 

REAtK  15  )(  COBE<  I  >, 1=1, NUMNP  ) 

REAIK  15  )(  XR(  I  )» 1=1 , NUMNP  ) 

REAIK  15  X  XT<  I  ),  1  =  1, NUMNP) 

REAIK  15  X  XZ<  I  ),  1=1,  NUMNP) 

REAIK  15  X  R(  I  ), 1  =  1, NUMNP  ) 

REAIK  15  X  Z<  I  ),  1=1,  NUMNP) 

REAIK  15  X  <  IX<  I*J)»J=1,5)»I  =  1*  NUMEL  ) 

REAIK  15  X  BETA<  I  ),  1  =  1, NUMEL) 

REAIK  15  X  ALPHA<  I  ),  1=1, NUMEL) 

REAIK  15  X  TEM(  I  ),  1=1 , NUMEL  ) 

DO  200  1=1, NUMMAT 

REAIK  15  )MTYPE»NT  »R0<  MTYPE  ) 

REAIK  15  X  <  E(  II*  J,  MTYPE),  J=l,  14),  11=1,  NT) 

DO  200  K=NT ,12 
DO  200  L=1 ,6 

200  E<  K,L, MTYPE  )=E<  NT,  L,  MTYPE) 

RETURN 

END 

SUBROUTINE  INTER 

COMMON/ ARG/RRRK  5  ),ZZZ<  5  ),RR<  4  ),ZZ<  4  ),S(  15,15  ),P(  15  >»TT(  6  ), 

1H(  6,15  ),CRZ<  6,6  ),XI<  10  ),  ANGLE<  4  ),SIG(  18),EP3<  18),N 
COMMON/PLANE/NPP 
DIMENSION  XM<  7)»R<7),Z(7)»  XX<  9  ) 

DATA  XX/3*. 1259391805448,3* . 1323941527884, .225, 

1  .696140478028, ,410426192314/ 

R<  7  )=(  RR(  1  HRR(  2  HRF<(  3  )  )/3 , 0 
Z(  7  )=(  ZZ(  1  )+ZZ(  2  )+ZZ<  3  )  )/3 . 0 
DO  100  1=1,3 
J=I+3 

R(  I  )=XX(  8  )*RR<  I  )+<  1 . 00-XX<  8  )  )*R<  7  ) 

R<  J  )=XX<  9  >*RR<  I  )+<  1 . 00-XX<  9  )  >*R<  7  ) 

Z(  I  )=XX(  8  )*ZZ<  I  )+<  1 ,00-XX<  8  )  >*Z<  7  > 

100  Z<  J  )=XX(  9  )*ZZ<  I  )+(  1  .OO-XXI  9  )  )*Z<  7  ) 

DO  200  1=1*7 
200  XM<  I  >=XX<  I  >*R<  I  ) 

DO  300  1=1,10 
300  XI<  I  >=0.00 

AREA= .  50*<  RR(  1  >*<  ZZ(  2  )-ZZ<  3  )  )+RR<  2  )*(  ZZ(  3  )-ZZ<  1  ) )  M<R(  3  )*<  ZZ(  1  ) 

1  -ZZ<  2  ))> 

IF( NPP.NE.O  )  GO  TO  600 

DO  400  1=1*7 

XI<  1  )=XI<  1  )+XM(  I  ) 

XI(  2)=XI(2)+XM(  I  )/R<  I  ) 
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XIC  3  >=XI<  3  HXMC  I  )/<  RC  I  H*2  ) 

XI(  4  )=XI(  4  HXMC  I  >*ZC  I  )/RC  I  ) 

XIC  5  )=XI(  5  HXMC  I  >*ZC  I  )/C  RC  I  )**2  ) 

XI <  6  >=XI(  6  HXMC  I  )*<  ZC  I  )#*2  >/<  F<«  I  )«2  ) 

XI<7)=XIC7HXMCI)*RCI> 

XIC8)=XIC3HXMC  I  >*ZC  I  ) 

XI(  9  >=XIC  9  HXMC  I  )*<  R<  I  )**2  ) 

400  XI<  10  )=XIC  10  HXMC  I  )*RC  I  )*Z<  I  ) 

no  500  1=1,10 

500  XI<  I  >=XI<  I  )*AREA 
RETURN 

600  XIC  1  )=AREA 

XI<  7  )=R<  7  >*AREA 
XIC  8  )=Z C  7  )*AREA 
RETURN 
END 

SUBROUTINE  MESH 
INTEGER  CODE 

DIMENSION  ARC  4,  8),AZ<  4,  8), NCODEC  4,  8) 

COMMON/TD/IMINC 100  ),IMAXC 100HJMINC 25 ) ,  JMAXC  25  > , MAXI » MAX J , NMTL , NBC 
COMMON/NPDATA/RC  10  ),  CODEC  10  >,XR(10  ),Z(10  ),XZ<10  )» 

1NPNUMC  4,  8  ),T<  10  ),XT<10  ) 

COMMON/ELDATA/BETAC  10  >,EPR<  10  >,PR(  4  ),SHC4  )»IX<8  ,5), 

1 1  PC  4  >,JPC4  ),  ISC  4  ),  JSC  4  ),ALPHAC  10  ),ITC4  ),JTC4  ), 

2STC  4  ) 

EQUIVALENCE  C  RC  1  )»AR)»C  ZC  1  ),AZ  ),C  IXC 1 , 1 ) , NCQDE  ) 


c* 

* 

************* 

*  *  *  * 

*  *  *  * 

* 

$ 

* 

* 

* 

* 

4 

* 

* 

4 

* 

* 

4 

c 

MESH  CONTROL  INFORMATION 

c* 

* 

************* 

*  *  *  * 

*  *  *  * 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

READ  C 5,1000)  MAXI,MAXJ, 

NSEG, NBC, 

NMTL 

URITEC 6,2000)  MAXI,MAXJ, 

NSEG, NBC, 

NMTL 

c* 

* 

************* 

*  *  *  * 

*  *  *  * 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

c 

INITIALIZE 

c* 

* 

^  it*  *  t*  ,«■  ^  iL  ib  ^  ^  1 1'  (X1  il' 

*  *  *  * 

*  *  *  * 

* 

n' 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

4 

ISEG=-1 
PI=3. 1415927 
DO  110  J=1 » 8 
DO  100  1=1,4 
NCODEC  I ,  J  )=0 
ARC  I,  J  )=0. 

AZC  I,  J  >=0. 

JMAXC  I  >=0 
100  JMIffC  I  )=MAXI 
IMINC  J  >=MAXJ 
110  IMAXC  J  )=0 

C#  ************************************ 
C  LINE  SEGMENT  CARDS 

C#  *******##*#*#*#***************&**** 
150  ISEG=ISEG+1 

159  IFC ISEG , EQ . NSEG  )  GO  TO  400 

READC  5,1001  )  I1,J1»R1»Z1,I2»J2»R2»Z2»I3»J3»R3»Z3»  IPTIQN 

WRITEC  6,2001  )I1 ,  J1  ,R1  ,Z1 , 12,  J2, R2, Z2, 13,  J3,R3, Z3,  IPTION 

IPTI0N=IPTI0N+1 

ARC  I1,J1  )=R1 

AZC  II, J1  )=Z1 

NCODEC I1,J1 )=1 

CALL  MNIMXC II , J1 > 

GO  TO  C 150,200,200,300,300,200,200 ),  IPTION 
C#  ***************************  ******** 
C  .  GENERATE  STRAIGHT  LINES  ON  BOUNDARY 
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L*  #*.«*****#«$$**  *********  ******* 

200  DI  =  ABS< FL0AT< 12-11  ) ) 

D.J=  ABS<  FL0AT<  J2-J1  )  ) 

AR<  I2»  J2  )=R2 
AZ<  1 2  *  J  2  )-Z2 
NC0DE(  12*  J2  )=1 
CALL  MNIMXC I2*J2  ) 

ISTRT=I1 

ISTP=I2 

JSTRT=J1 

JSTF-J2 

DIFF=MAX1<  DI*BJ  ) 

ITER=DIFF-1 ♦ 

I INC=0 
JINC=0 

IF(  I2.NE.il  )  IINC=< 12-11  )/IABS(  12-11) 

IF(  J2.NE.J1 >  JINC=(  J2-J1  )/IABS< J2-J1 ) 

KAPPA=1 

I F(  1 2 . ME . 1 1 . AND . J2.NE.J1 . AND , IPTI ON . NE . 3  )  KAPPA=2 
IF(  KAPPA. EQ.2  )  DIFF=2 . *DIFF 
RINC=(  R2-R1  )/DIFF 
ZINC=( Z2-Z1  )/DIFF 

WRITE(  6*2002  )  DI* DJ»DIFF * R INC* ZINC* ITER* I INC* JINC* KAPPA 
C 

C  CHECK  FOR  INPUT  ERROR 

C 

IF< KAPPA.NE.2.0R.DI.EQ.DJ  )  GO  TO  210 
URITE<  6*2003  ) 

GO  TO  150 
C 

C  INTERPOLATE 

C 

210  1=11 
J=J1 

WRITE( 6*2004) 

DO  230  M=1 » ITER 

IF( ITER. EQ.O. AND. IPTION. EQ.2)  GO  TO  230 

IF(  ITER • EQ . 0. AND • IPTI ON .EG . 6  )  GO  TO  230 

IF< ITER. EQ.O. AND. IPTION. EQ. 7  )  GO  TO  230 

IF(  KAPPA. EQ.2)  GO  TO  220 

IOLD=I 

1=1+1 INC 

JOLD=J 

J=J+.JINC 

AR(  I  *  J  )=AR<  IOLD* JOLD  )+RINC 
AZ(  I  * J  )=AZ(  IOLD* JQLB  )+ZINC 
URITE< 6*2005)  I  *  J»AR<  I » J  )*  AZ<  I  *  J  ) 

CALL  HNIMX< I»J  ) 

NCODE<  I»J)=1 
GO  TO  230 
220  CONTINUE 

IF< II. GT. 12. AND. IPTION. EQ. 7  )  GO  TO  221 
IF< II. LT. 12. AND. IPTION. EQ. 6)  GO  TO  221 
IOLD-I 
1=1+1 INC 

AR< I »J  )=AR< IOLD* J  HRINC 
AZ<  I » J  )=AZ<  IOLD*  J  )+ZINC 
WRITE<  6*2005  )  I*J*AF«  I*J)*AZ<  I.J) 

NCODE<  I  *  J  )=1 
CALL  MNIMX( I*J  ) 
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/ 


non 


i 


JOLD=J 
J= J+ JING 

AR<  I  »J  )~AR(  If JOLD  HRINC 
AZ<  I » J  )=AZ<  I » JOLD  HZINC 
NCODEC If  JH1 

WRITE* 6»  2005  )  I » J » AR( I » J  ) >  AZ< I » J ) 

CALL  NNIMX* I » J  ) 

GO  TO  230 
221  J0LD=J 
J=JhJINC 

AR(  If  J)=AR<  If  JOLD  HRINC 
AZ(  If  J)=AZ<  If  JOLD  HZINC 
NC0DE(  If  J  )=1 

WRITE*  <S»  2005  )  If.Jf  AR(  If  J  )fAZ<  If  J  > 

CALL  MNIMX(IfJ) 

I0LD=I 

I=I+IINC 

AR(  If  J)=AR<  IOLDf  J  HRINC 
AZ<  If  J  )=AZ<  IOLDf  J  HZINC 
NCODE* If J  )=1 

WRITE( 6f  2005  )  If J»AR( If J  )fAZ* I»J  ) 

CALL  HNIMX* If J  ) 

230  CONTINUE 

IF< KAPPA. EQ. 1  )  GO  TO  150 

IF(  II  »GT .  I2.ANEI.  IPTI0N.EQ.7  )  GO  TO  231 

IF(  I1.LT.I2.AND.IPTI0N.EQ.6)  GO  TO  231 

IOLD=I 

I=I+IINC 

AR<  I  f  J  )=AR*  IOLDf  J  HRINC 
AZ(  If  J  )=AZ(  IOLDf J  HZINC 
GO  TO  232 

231  CONTINUE 
JOLD=J 
J=J+.JINC 

AR<  I  f  J  )=AR<  If  JOLD  HRINC 
AZ(  IfJ)=AZ<  If  JOLD  HZINC 

232  CONTINUE 
NCODEC  If  J  >=1 

WRITE*  6f  2005  )  I  f  J,  AR<  I  f  J  )f  AZ<  If  J  ) 

CALL  MNINX(IfJ) 

GO  TO  150 

C*  *********************************** 
C  GENERATE  CIRCULAR  ARCS  ON  EOUNDARY 

C*  *********************************** 
300  AR<  I2fJ2)=R2 
AZ(  1 2  f  J 2  )=Z2 
NCODE( I2f J2 )  =  1 
CALL  MNIMX<I2fJ2) 

IF( IPTI0N.EQ.5  )  GO  TO  320 

FIND  CENTER  OF  CIRCLE 

AR<  I3f  J3  )=R3 
AZ<  I3f  J3  )=Z3 
NCODE( I3f J3 )=1 
CALL  HNIMX<I3fJ3) 

SLAC=( Z2-Z1 >/(  R2-R1 ) 

SLBF=-1,/SLAC 
SLCE=<  Z3-Z2 )/<  R3-R2 ) 

SLDF=-1 . /SLCE 
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c 

C  CHECK  FOR  INPUT  ERROR 
C 

IF<  ABS(  SLAC-SLCE  ).GT  . .001  )  GO  TO  310 
UJRITE<  6.2006)  R1.Z1.R2.Z2.R3.Z3.SLAC.SLCE 
GO  TO  150 

310  R4=Rl+( R2-R1  )/2. 

Z4=Z  1  +(  Z2-Z1  )/2. 

R5=R2+<  R3-R2  )/2. 

Z5=Z2+(  Z3-Z2V2. 

BBF=Z4-SLBF*R4 
BDF=Z5-SLDF*R5 
RC=(  BBF-BDF  )/<  SLBF-SLBF  ) 

ZC=SLBF*RC+BBF 
WRITE< 6.2007)  RC.ZC 
KAPF’A=1 
GO  TO  330 
320  KAPPA=2 
RC=R3 
ZC=Z3 

330  ISTRT=I 1 
ISTP=I2 
JSTRT= J 1 
JSTF'=J2 
RSTRT=R1 
RSTP=R2 
ZSTRT=Z1 
ZSTP=Z2 

340  CALL  ANGLE<  RSTRT .ZSTRT . RC.ZC. ANG1 > 

CALL  ANGLE< RSTP , ZSTP .RC.ZC.ANG2 ) 

IF( ANG2.LE. ANG1  )  ANG2=2.0*PI+ANG2 
C 

C  FIND  ANGULAR  INCREMENT 
C 

DI=  ABS< FLOAT< ISTP-ISTRT  )> 

D.J=  ABS(  FLOAT<  JSTP-JSTRT  > ) 

I INC=0 
JINC=0 

IF< ISTRT.NE.ISTP  )  IINC=( ISTP-ISTRT  >/IABS< ISTP-ISTRT) 
IF( JSTRT.NE. JSTP  )  JINC=<  JSTP-JSTRT  )/IABS(  JSTP-JSTRT ) 
LAMDA=1 

IF(  I INC . NE <0. AND . JINC  *NE , 0  )  LAMBA=2 
0IFF=MAX1< DI.DJ  ) 

ITER=DIFF-1 ♦ 

IF( LAMDA.EQ.2)  BIFF=2.*BIFF 
DELPHI=( ANG2-ANG1  )/DIFF 
WRITE( 6.2008  )  ANG1 . ANG2 . BIFF .DELPHI 
C 

C  CHECK  FOR  INPUT  ERROR 

C 

IF<  LAMDA»NE.2*0R*BI tEQ.DJ  )  GO  TO  350 
WRITE< 6.2003  ) 

GO  TO  150 
350  I0=ISTRT 
JO=JSTRT 
WRITE<  6.2004) 

C 

C  INTERPOLATE 

C 

NPT=IABS<  12-11  )+IA8S(  J2-J1  )-l 
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LIU  330  N=1  -i  1 ILK 
35?  IF!  L.AMDA .  EQ . 2 )  GO  TO  360 
'  1  =  10  PI INC 

J=JO+JINC 
CALL  MN1NX! I » J  > 

NCODE!  I » J  )=1 

CALL  CIRCLE< ANG1 » DELPHI * RSTRT » ZSTR T # RC » 1C , I # J > 

WRITE!  6  y  2005  )  I  *  J»AR!  I*  J  >»AZ<  I,  J  ) 

GO  TO  370 
360  I=IO+IINC 
J=JO 

NCODE<  I»  J  >=1 
CALL  MNIMX! I  *  J  ) 

CALL  CIRCLE! ANG1 » DELPHI * RSTRT * 2STRT * RC» 2C *  I » J  ) 

WRITE! 6  * 2005  )  I » J  *  AR! I » J  ) » AZ( I » J  ) 

J=JO+JINC 
NCODE!  I » J  )=1 
CALL  MNIhX! I*J  ) 

CALL  CIRCLE! ANG1 , DELPHI * RSTRT * ZSTRT .  RC . ZC . I >  J  > 

WRITE!  6,2005)  I,J»AR!  I*J  >*AZ!  I»J> 

370  10=1 
330  JO=J 

IF! LAMDA.NE.2  )  GO  TO  390 
I=IO+IINC 
NCODE!  I  *  J  )=1 
CALL  MNIHX! I » J  ) 

CALL  CIRCLE!  ANG1 , DELPHI , RSTRT * ZSTRT *  RC  *  ZC , I , J  > 

WRITE!  6*2005)  I,J*AR!  I*J  )»AZ!  I,.J  ) 

390  IF! KAPPA. EG. 2)  GO  TO  150 
ISTRT=I2 
I STP=I3 
JSTRT=J2 
JSTP=J3 
RSTRT=R2 
RSTP=R3 
ZSTRT=Z2 
ZSTP=Z3 
KAPPA=2 

399  GO  TO  340 

r<if  ^  UL-  ^  \L-  d 

wT  ^  ^  *5^  ^  ^  *  *T  x  *1^ 

C  CALCULATE  COORDINATES  OF  INTERIOR  POINTS 

C#  #***###***#****###*  *$***********&** 

400  IF!  HAXJ  *LE . 2  )  CO  TO  430 
J2=MAX J-l 

DO  420  N=1 » 500 
RESID=0. 

DO  410  J=2 ,  J2 
1 1  =  IMIN!  J  HI 
I2=IHAX! J  )-l 
DO  410  1=11*12 

IF!  NCODE! I  * J  )*EQ. 1  )  GO  TO  410 

DR=!  AR!  IP1*J  )+AR!  I~1*J  HAR!  I,  J+l  HAR!  I  *  J -1 )  )/4 . -AR!  I » J ) 

DZ=!  AZ!  IP1  *  J  )+AZ!  I-1,J  HAZ!  I»JH  HAZ!  I*  J-l  ) )/4 ,-AZ< I* J ) 
RESID=RESID f  ABS!  DR  HABS!  DZ  ) 

AR!  I,J )=AR!  I»J  HI « 8#DR 
AZ!  I*J  )=AZ!  I  *  J  H1.8*DZ 
410  CONTINUE 

IF!  N»EQ« 1  )  RES1=RESID 

IF!  N*EQ» 1 ♦ AND »RESID *ECK 0 ♦  )G0  TO  430 

IF! RESID/RES1 »LT  *  1 *E-5  )  GO  TO  430 
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420  CONTINUE 

430  WRITE( 6 r 2009 )  N 

WRITE< 15)NBC»NMTL 

“C*  *****************  ****************** 
CALL  POINTS 

•  C*  *********************************** 
1000  FORMAT  (515) 

100.1.  FORMAT  < 3< 2I3.2FG.3 )»I5 ) 

2000  FORMAT  < 30H1  MESH  GENERATION  INFORMATION// 


1  41  HO  MAXIMUM  VALUE  OF  I  IN  THE  MESH - 13/ 

2  41  HO  MAXIMUM  VALUE  OF  J  IN  THE  MESH - 13/ 

3  41H0  NUMBER  OF  LINE  SEGMENT  CARDS - 13/ 

4  41H0  NUMBER  OF  BOUNDARY  CONDITION  CARDS - 13/ 

5  41H0  NUMBER  OF  MATERIAL  BLOCK  CARDS - 13/// > 


2001  FORMAT  ( //8SH  INPUT  II  J1  R1  Z1  12  J2  R2  Z 

12  13  J3  R3  Z3  IPTI0N/8X»3(2I4,2F8.4  )rI6 ) 

2002  FORMAT  ( 5H  DI=F4.0»5H  DJ=F4.0,7H  DIFF=F4.0»7H  RINC=F8.3»7H  ZI 

1NC=F8.3»7H  ITER=I3»  7H  IINC=I3»7H  JINC=I3»8H  KAPPA=I1) 

2003  FORMAT<  1X»38H**BAD  INPUT— THIS  LINE  IS  NOT  DIAGONAL) 

2004  FORMAT  < 30H  I  J  AR  AZ ) 

2005  FORMAT  <2I5»2F11.6) 

2006  FORMAT  < 51H  **  BAD  INPUT  -  THESE  POINTS  DO  NOT  DEFINE  A  CIRCLE»/» 
13Xr6F12.4» 10X»2E20.8  ) 

2007  FORMAT(  19H  CENTER  COORDINATE ,< Fll .6» 1X»F11 . 6»1X  ) ) 

2008  FORMAT  ( 7H  ANG1=F9.6.7H  ANG2=F9.6»7H  DIFF=F3.0»9H  DELPHI=F9.6) 

2009  FORMAT  < //30H  COORDINATES  CALCULATED  AFTER  I3»11H  ITERATIONS) 
RETURN 

END 


SUBROUTINE  MINV 
PURPOSE 

INVERT  A  MATRIX 
USAGE 

CALL  MINV<  Af Nf D»L»M  ) 

DESCRIPTION  OF  PARAMETERS 

A  -  INPUT  MATRIX?  DESTROYED  IN  COMPUTATION  AND  REPLACED  BY 
RESULTANT  INVERSE. 

N  -  ORDER  OF  MATRIX  A 
D  -  RESULTANT  DETERMINANT 
L  -  WORK  VECTOR  OF  LENGTH  N 
M  -  WORK  VECTOR  OF  LENGTH  N 

REMARKS 

MATRIX  A  MUST  BE  A  GENERAL  MATRIX 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 

METHOD 

THE  STANDARD  GAUSS-JORDAN  METHOD  IS  USED.  THE  DETERMINANT 
IS  ALSO  CALCULATED.  A  DETERMINANT  OF  ZERO  INDICATES  THAT 
THE  MATRIX  IS  SINGULAR. 


SUBROUTINE  MINV<  A.N»D»L»M  ) 
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DIMENSION  A(  1  )»L<  1  )  »M(  1  ) 


IF  A  DOUBLE  PRECISION  VERSION  OF  THIS  ROUTINE  IS  DZ3IKEB,  THE 
C  IN  COLUMN  1  SHOULD  BE  REMOVED  FROM  THE  DOUBLE  PRECISION 
STATEMENT  WHICH  FOLLOWS. 

DOUBLE  PRECISION  A,D»BIGA,HOLD 

THE  C  MUST  ALSO  BE  REMOVED  FROM  DOUBLE  PRECISION  STATEMENTS 
APPEARING  IN  OTHER  ROUTINES  USED  IN  CONJUNCTION  WITH  THIS 
ROUTINE. 

THE  DOUBLE  PRECISION  VERSION  OF  THIS  SUBROUTINE  MUST  ALSO 
CONTAIN  DOUBLE  PRECISION  FORTRAN  FUNCTIONS.  ABS  IN  STATEMENT 
10  MUST  BE  CHANGED  TO  DABS. 


SEARCH  FOR  LARGEST  ELEMENT 

D=1 .0 
NK=-N 

DO  80  K=1  ,N 
NK=NK+N 
L(  K >=K 
M(  K  )=K 
KK=NK+K 
BIGA=A( KK  ) 

DO  20  J=K»N 
IZ=N*< J-l ) 

DO  20  I=K*N 
IJ=IZ+I 

10  IF(  ABS(  BIGA )-ABS< A(  IJ  ) )  )  15,20,20 
15  BIGA=A<  I J  ) 

L(K  )=I 
M(  K  >=J 
20  CONTINUE 

INTERCHANGE  ROWS 


J=L<  K  ) 

IF(J-K)  35,35,25 
25  KI=K-N 

DO  30  1  =  1, N 
KI=KI+N 
HOLD=-A(  KI ) 

JI=KI-K+J 
A(  KI  )=A(  JI  ) 

30  A< JI  )  =HOLD 

INTERCHANGE  COLUMNS 

35  I=M<  K  ) 

IF(I-K)  45,45,38 
38  JP=N#< 1-1  ) 

DO  40  J=1,N 

JK=NK+J 

JI=JP+J 
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HOLD=-A<  JK ) 

A<  JK  )=A(  JI  > 

>  40  A<  JI  )  =H0LD 

DIVIDE  COLUMN  BY  MINUS  PIVOT  (VALUE  OF  PIVOT  ELEMENT  IS 
CONTAINED  IN  BIGA ) 

45  IF< BIGA  )  40 >46*48 

46  D=0.0 
RETURN 

48  DO  55  1=1  »N 

IF(I-K)  50 » 55 » 50 
50  IK=NK+I 

A(  IK  )=A<  IK  )/<  -BIGA  ) 

55  CONTINUE 

REDUCE  MATRIX 

DO  65  1=1  »N 
IK=NK+I 
HOLB=A<  IK  ) 

I J=I-N 
DO  65  J=1  »N 
I J=I J+N 

IF(  I-K)  60 » 65 r 60 
60  IF(J-K)  62 » 65 » 62 
62  KJ=I J-I+K 

A<  I J  )=H0LD*A<  KJ  )+A<  I J  > 

65  CONTINUE 

DIVIDE  ROW  BY  PIVOT 

KJ=K-N 
DO  75  J=1 »N 
KJ=KJ+N 

IF(J-K)  70»75»70 
70  A(  KJ  )=A(  KJ  )/BIGA 
75  CONTINUE 

PRODUCT  OF  PIVOTS 

D=D*BIGA 

REPLACE  PIVOT  BY  RECIPROCAL 

A< KK  )=1.0/BIGA 
80  CONTINUE 

FINAL  ROW  AND  COLUMN  INTERCHANGE 

K=N 

100  K=<  K-l  ) 

IF<  K  )  150» 150» 105 
105  I=L<  K ) 

IF< I-K  >  120f 120» 108 
108  JQ=N*(K-1  ) 

JR=N*( 1-1  ) 

DO  110  J=1 »N 
JK=JQ+J 
HOLD=A( JK  ) 
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A!  JK  >=-A!  JI  ) 

110  A!JI>  -HOLD 
120  J=M!  K  ) 

IF!  J~K  >  1 00 . 100, 125 
125  KI“K“N 

DO  130  I - 1  *  N 
KI=KI  FN 
HGLD=A!  KI  ) 

JI=KI-KiJ 
A!  KI  )=-A! JI  ) 

130  A! JI  )  -HOLD 
GO  TO  100 
150  RETURN 
END 

SUBROUT I NE  MN I MX!  I  *  J  > 

CQMMON/TD/IMIN< 100  ).IMAX! 100  ),JMIN!  25  ). JMAX! 25  ) *  MAXI » MAX J » NMTL  » NBC 

IF!  J.LT .  JMIN!  I  ) )  JMIN!  I  )=J 

IF(  J.GT.JMAX! I >)  JMAX! I  )=J 

IF<  I.LT.IMIN!  J))  IHIN(  J  )  =  I 

IF<  I  »GT  .  IMAX!  J  ) )  IMAX(J)=I 

RETURN 

END 

SUBROUTINE  MODIFY!  NEQ.N.U  ) 

COMMON/SOLVE/B!  72).  A!  72. 36 > > NUMBLK, MBAND 

DO  10  M=2» MBAND 

K=N-MT1 

IF<  K.LE.O  )  GO  TO  5 
B!  K  )=B!  K  )-A(  K.M  )*U 
A!  K.M  )=0.00 
5  K=N+M-1 

IF<  NEQ.LT .K  )  GO  TO  10 
B(  K  )=B!  K  )-A(  N .  M  >«U 
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A(  N.M  >=0.00 
CONTINUE 
A!  N»  1  >=1.00 
B(  N)=U 
RETURN 
END 

SUBROUTINE  MPLOT 
INTEGER  CODE 

COMMON/TD/IMIN!  100  >. IMAX< 100  ). JMIN!  25  >, JMAX! 25  ) . MAXI . MAX J » NMT1 
COMMON/NPDAT  A/R!  10  >.  CODE!  10  )?XR!  10  )»Z!  10  >.XZ<  10  >. 

1NPNUM!  4.  8  )  >  T!  10  ).XT!10  > 

REAL  X!  100  ).Y!  100  )»TX!  2  )»TY!  2  )» TITLE!  20  ).  ZMAX 
!  5.1000)  TITLE.RMAX.ZMAX 
CCP2SY  ! 0.7 .0.2. 0,2. TITLE. 0.0. 30 > 

CCP1PL  !  0,7.0. 7. --3  )  ^ 


READ 

CALL 

CALL 


TX!  1  )=0,00 


TY!  1  >=0,00 

TX!  2  )=RMAX/9.0 

TY!  2  )=RMAX/9»0 

ZMAX=ZMAX#TY!2H2.0 

IF  ! ZMAX.LT .17.0)  ZMAX=17.0 


'DO  100  J=1 »MAXJ 
NSTART=IMIN!  J ) 
NSTOP=IMAX!  J) 

N=0 

DO  101  I=NSTART .NSTOP 
N=N+1 
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NP=NPNUM<  I,.J  ) 

Y<  N  )=R<  NP  ) 

'  101  X<  N  )=Z(  NP  > 

C  CALL  CCP6LN  < X*Y,N* 1 *TX,TY ) 

100  CONTINUE 

DO  102  1  =  1  *  MAXI 

NSTART=JMIN< I  )  • 

NSTOP=JMAX< I  ) 

N=0 

DO  103  J=NSTART,NSTQP 
N=N+1 

NF‘=NPNUM<  I, J  ) 

Y<  N  >=R<  NP  ) 

103  X<  N  )=Z<  NP  ) 

C  CALL  CCP6LN  <  X* Y * N* 1 *TX , TY ) 

102  CONTINUE 

C  CALL  CCPIF'L  < ZMAX , -0 . 7 , -3  ) 

1000  FORMAT  ( 20A4/2F10 . 0  ) 

RETURN 

END 

SUBROUTINE  NAXSTF< 1 1 » J J*  KK  > 

INTEGER  CODE 

COMMON/ VISC/EF‘SDN<  12* 10  *  8  )  * SIGVPv  6  >.DEPSR<  6*10.8  )» BELT  IN 
COMMON/PLAS/ALFA<  6.  4 * 8  )*SIGYLDC  7*6*8)*  IFGF'L(  4*8) 

COMMON/NXDATA/NTP » NTS , NTOTS. GTS1G(  24.24,8) 

COMMON/MATP/RO<  6  ).E(  12.16,6  ),EE(  16  )*A0FTS(6) 

COMMON/BASIC/ACELZ ,  ANGVEL ,  ANGACC ,  TREF  *  VOL ,  NUMNP ,  NUMEL.NUMPCtNUMSC , 
1NUMST 

COMMON/ ARG/RRR(  5  ) ,  ZZZ<  5  > » RR<  4  ) ,  ZZ<  4  )  *  S(  15 , 15  ) »  P(  15  ) ,  TT(  6 )  * 

1H(6*15),  CRZ(  6*6), XI(  10)*ANGLE(  4  ),SIG(  18)*EPS<  18),  N 
COMMON/NPDATA/R<  10  ),C0DE(10  ),XR(10  )*Z(10  ),XZC10  ’)* 

1NPNUM(  4,  8 ), T<  10  )*XT(  10  ) 

CGMMON/ELDATA/B£TA(  10  ),EF*R<  10  )»PR<4  ),SH(4  )*IX(8  *5), 

1  IP(  4  ),  JP<  4  ),  IS(  4  )  *  JS<  4  ),  ALF’HA(  10  ),IT<4  )*JT(4  ), 

2ST<  4  ) 

COMMON/NXQUAD/AR 1 

COMMON/NONAXI/SK  30*30  )*P1( 30  ) .THETA, B31C  6,30 ) 

DIMENSION  C<  13,18  ),B(  18,18),B1(  6,18),B2<  6, 13  ),B3C  6*  13  ),B4(  6 *  18  )* 

1  B5<  6 , 1 8  )  *  B6<  6  *  1 8  )  *  B 1  A(  6 , 1 3  )  *  B1BC  6 , 18  ) ,  B2AC  6 , 1 8  )  *  B2B(  6  *  18 

2  )*B3A(  6*  18  )*B3B(  6,18  ),B4A(  6*18  ) ,B4B(  6*  18  ) *B5A(  6*18), 

3  B5B<  6*18  )*B6A( 6*18  )*B6B<  6,18  )*TVF( 18  ) 

C  ZERO  MATRICES 

DO  100  1=1*13 

DO  100  J=1 *18 

100  C<  I,J  )=  0*0 
DO  110  1=1*6 
DO  110  J=1 *  18 
Bl(  I*J)  =0.0 
B2( I*J )  =0.0 
B3<  I  *.J  >  =0,0 
B4(  I,J)  =0.0 
B5( I,J )  =0.0 
110  B6< I*J  )  =0.0 

RR<  1  )  =  RRR(  II  ) 

RR(2)  =  RRR(JJ) 

RR<  3  )  =  RRR(KK) 

ZZ(  1  )  =  ZZZ<  II  ) 

ZZ(  2  )  =  ZZZ(JJ) 

ZZ<  3  )  =  ZZZ(KK) 

COMM=RR<  2  )*<  ZZ<  3  )-ZZ<  1  )  )+RR(  1  )*<  ZZ(  2  )-ZZ(  3  )  )  FRR(  3  )U  ZZ(  1  )-ZZ(  2  )  ) 
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C  FILL  C  INVERSE 

C< 1 x  1  )=  (  RR<  2  )*ZZ<  3  )  -RR<  3  )*  ZZ<  2  ) )  /  COMM 
„  C<  1 x4  )=  (  RR(  3  )*ZZ<  1  )  -RR(  1 )»  ZZ<  3  ) )  /  COMM 

C<lx7)=  <  RR(1)*ZZ<2)  -RR<2)*  ZZ(  1  ) )  /  COMh 
C<  2x  1  )=  <  ZZ<  2  )  -  ZZ<  3  )  >  /  COMM 
-  C<  2* 4  )=  <  ZZ<  3  )  -  ZZ(  1  ) )  /  COMM 

C(  2  x  7  )=  <  ZZ<  1  )  -  ZZ<2))  /  COMM 

C<  3x  1  )=  <  RR<  3  )  -  RR<  2  )  )  /  COMM 

C<  3x 4  )=  <  RR(  1  )  -  RR(  3  ) )  /  COMM 

C<3x7)=  <  RR<  2  )  -  RR(  1  ) )  /  COMM 

C<  4x2  )=  C(  1.1  ) 

C<  4x5  )=  C(l»4) 

C<  4  x  8  >=  C<1*7) 

C<  5»2  )=  C<  2»  1  ) 

C<5x5)=  C<2x4) 

C<  5x8  )=  C(  2x7  > 

C<6x2)=  C(  3 x  1  ) 

C<  6x5  >=  C<  3x4  ) 

C(  6x8  )=  C(  3x7  ) 

C(  7x3)=  C(  lxl  ) 

C(  7x6  )=  C<  1x4) 

C(7x?)=  C(  1x7) 

C<  8x3  )=  C(  2x  1  ) 

C(  8x6  )=  C<  2x4  ) 

C<  8x 9  )=  C<  2x7  ) 

C(  9x3  )=  C<  3x1  ) 

C<  9x6  )=  C<  3x4  ) 

C(  9x9  )  =  C( 3x7  ) 

DO  120  1=10x18 
DO  120  J=  1x9 
II  =  1-9 
Jl=J+9 

C<IxJ)  =(-l. /THETA)  *  C<IlxJ) 

C(  I  x  J1  )=<  1.  /THETA)  *  C<IlxJ) 

120  CONTINUE 
C  FILL  B  MATRICES 
C  B1  CONSTANT  TERMS 
C  B2  THETA  TERMS 
C  B3  1/R  TERMS 
C  B4  THETA/  R  TERMS 
C  B5  Z/R  TERMS 
C  B6  THETA  *Z/R  TERMS 
DO  130  J=1 x 18 
Bl<  1  x  J  )  =  C(2xJ) 

Bl(  2x  J  )  =  C<6xJ) 

Bl<  3x  J  )  =  C(2xJ)+C(  17 x  J  ) 

Bl<  4x  J  )  =  C(  3x  J  )+C<  5xJ  ) 

Bl(  5x  J  )  =  C(9xJ)  +C<  14x  J  ) 

Bl(  6x  J  )  =  C(llxJ) 

B2<  lxJ)  =  C(llxJ) 

B2(  2x  J  )  =  C(15xJ) 

B2(  3x  J  )  =  C(llxJ) 

B2<  4x  J  )  =  C<  12xJ  )+C(  14x  J  ) 

B2<  5x  J  )  =  C<  18x  J  ) 

B3<  3x  J  )  =  C(  lxJ  )+  C<  1 6 x  J  ) 

B3<  5x  J  )  =  C<13xJ) 

B3(  6x  J  )  =  C(10xJ)  -  C<  7 x  J  ) 

B4<  3x J  )  =  C< lOx J  ) 

B4<  6x  J  )  =  -C(  16x  J  ) 

B5(  3 x  J  )  =  C(  3 x  J  )  +C(  18x  J  ) 
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B5<  5 » J  )  =  C(  15.  J  ) 

B5<  6,  J  )  =  C<  12  >J  >~C(  9».J  ) 

B6(3,,J>  =  C<  1 2 » J  ) 

B6(  6,  J  )  =  -C<  18,  J  ) 

130  CONTINUE 

C  NOW  CALCULATE  BT  *  D  *  B 
CALL  INTER 

THETA2  =  < THETA  **2)/2.0 
THETA3  =  (THETA  **3)/3.0 
BO  140  1=1 » 6 
DO  140  J=1 » 18 

B1A<  I,J  )=<B1<  I,J  >*XI(  1  >  +B3(I,J>*  XI(  2  )  +  B5<I»J>*  XI(  4  )  >*  THETA  + 

1  (B2(I»J)*XI<1>  +B4<I,J>*  XI(  2  )  +  B6(I,J>*  XI(  4  )  )*  THETA2 

B2A(  I » J  >=(  B1(I»J)*XI<1)  +B3(  I ,  J  )*  XI(  2  )  +  B5(I,J>*  XI(  4  )  )*  THETA2 

1  +  <B2<  I,J)*XI<  1  )  +B4(I,J)*  XI(  2  )  +  B6(I»J>*  XI(  4  )  >*  THETA3 

B3A<  I,J)=(B1<  I,J)*XI<2)  +B3(  I » J  )*  XI(  3  )  +  B5(I,J>*  XI<  5  >  >*  THETA 

1  +<  B2<  I,J)*XI(2)  +B4(I,J>*  XI(  3 )  +  B6(I,J>*  XI<5>>*  THETA2 

B4A(  I  » J  )=<  Bl<  I  *  J  )*XI<  2  )  +B3(  I » J  )*  XI(  3  )  +  B5(I.J>*  XI(5>>*  THETA2 

1  +<  B2<  I,J)*XI<2)  +B4(I,J)*  XI(  3 )  +  B6(I,J)*  XI<5)>*  THETA3 

B5A<  I  »  J  )=(  Bl(  I » J  )*XI<  4  )  +B3(I,J>*  XI(  5  )  +  B5(I,J)*  XI(6>>*  THETA 

1  +  <  B2<  I,J  >*XI<4  )  +B4<I»J>*  XI(  5  )  +  B6(I,J>*  XI(  6  )  )*  THETA2 

B6A<I»J)=<B1<I,J>*XI<4>  +B3<  I » J  )*  XI(  5  )  +  B5(I»J>*  XI(6>>*  THETA2 
1  +  (  B2<  I,J  >*XI<4  >  +B4(I,J>*  XI(  5  >  +  B6(I,J>*  XI(  6  )  )*  THETA3 

140  CONTINUE 

DO  150  1=1 ,6 
DO  150  K=1 » 18 
B1B<  I  »K  )=  0.0 
B2B(  I  »K  )=  0.0 
B3B(  I  »K  )=  0.0 
B4B<  I,K  )=  0.0 
B5B(  I  »K  )=  0.0 
B6B<  I  »K  )=  0.0 
DO  150  J=l,6 

B1B(  I » K  )  =  B1B<  I  »K  )  +  CRZ(I»J)  *  B1A(J,K> 

B2B<  I  »K  )  =  B2B<  I  »K  )  +  CRZCIpJ)  *  B2A(.J,K> 

B3B<  I »K  )  =  B3B(I,K>  +  CRZ(I,J>  *  B3A(J,K> 

B4B<  I »K  )  =  B4B(  I » K  )  +  CRZ(I,J>  *  B4A(J,K> 

B5B<  I  »K  >  =  B5B(  I  »K  )  +  CRZCIpJ)  *  B5A(J,K) 

B6B(  I  »K  )  =  B6B(  I  ,N  )  +  CRZ(I,J>  *  B6A(J,K) 

150  CONTINUE 

DO  160  1=1,18 

DO  160  K=1 *  18 

B(  I ,  K  )=0,0 
DO  160  J=1 » 6 

B(I,K>  =  B<  I »K  )  +  Bl(J.I)*  B1B<  J,K)+B2(  J,I)*B2BC  J,K)+B3(J,I  )* 

1  B3B(  J»K  )+B4(  J,I  )*B4B<  J,K  )+B5(  J,  I  >*B5B(  J  ,  K  )+B6<  J » I  )*B6B<  J,K> 
160  CONTINUE 
250  CONTINUE 

C  B(ItK)  NOW  CONTAINS  THE  STIFFNESS  MATRIX  FOR  ONE  TRIANGULAR  ELEMENT 
AR1  =  AR1  +  XI( 1  )  *THETA 
DO  235  K=1 » 6 
DO  235  1=1,3 

BS1( K, 3*11-3+1  )  =  BS1(K, 3*11-3+1  )  +B1A(K,I  ) 

BS1<  K,3*.JJ-3+I  )  =  BSl(K,3*JJ-3+I  )  +BlA(K,I+3  ) 

BS1<  K » 3*KK-3+I  )  =  BS1< K,3*NK-3+I >  +BlA(K,I+6  ) 

BS1< K, 3*11+1+12)=  BS1(K, 3*11  +  12+1  )+BlA(K,I+?  ) 

BS1(  K,3*.J.J  +  I  +  12  >=  BS1<  K , 3*J J  +  12  +  I  )+BlA(  K,I+12) 

BS1(K,3*KK+I  +  12)=  BS1(K,3*NN+12+I  )+BlA(  N ,  1+15  ) 

235  CONTINUE 

IIM  =  3*  II  -3 
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JJM  =  3*  JJ  -3 
KKM  =  3:J<  KK  -3 
DO  170  K=1 » 4 
DO  170  1=1,3 
DO  170  J=1 » 3 

IF<  K.EQ. 1  .OR,  K.EQ.2)  11  =  1 
IF(K.EQ,3  .OR,  K.EQ.4  )  11  =  1  +? 

IF(  K.EQ.l  .OR.  K.EQ.3)  J1  =  J 
I F(  K.EQ.2  .OR.  K.EQ.4)  J1=.J  +9 

IF< K.EQ.l  .OR,  K.EQ.2)  K1=0 
IF(  K.EQ.3  .OR.  K.EQ.4)  Kl  =  15 
IF(  K.EQ.l  .OR.  K.EQ.3)  K2=0 
IF(  K ,EQ , 2  .OR.  K.EQ.4)  K2=15 
182  KK2=KKM 
II2=IIM 
JJ2=JJM 
180  KK1=KKM 
JJ1=JJM 
II1=II« 

Sl< II1+I+K1»II2+J+K2)  =  Sl< II1+I+K1*II2+J+K2)  +B<I1»J1> 

Sl< II1+I+K1,JJ2+J+K2>  =  Sl(  IIl+I+Klf  JJ2+J+K2)  +B(Il»Jl+3> 

Sl<  II1+I+K1,KK2+J+K2>  =  Sl< II1+I+K1 *KK2+J+K2 >  +B(Il,Jl+6> 

Sl< JJ1+I+K1,II2+J+K2>  =  Sl< JJ1+I+K1,II2+J+K2)  +B(I1+3,J1) 

Sl< JJ1+I+K1,JJ2+J+K2>  =  Sl< JJ1+I+K1,JJ2+J+K2>  +B< 11+3, Jl+3 > 

Sl<  JJ1+I+K1»KK2+J+K2)  =  Sl<  JJ1+I+K1 »KK2+J+K2  )  +B( 11+3, Jl+6 > 

Sl<  KK1+I+K1 » I I2+J+K2  )  =  Sl<  KK1  +  I+K1 » I I2+J+K2  )  +  B(I1+6,J1> 

Sl< KK1+I+K1,JJ2+J+K2)  =  Sl( KK1+I+K1 » J J2+ J+K2  )  +  B(Il+6,Jl+3> 

Sl<  KK1+I+K1 , KK2+J+K2  )  =  Sl< KK1+I+K1 , KK2+J+K2 >  +  B(Il+6,Jl+6> 

170  CONTINUE 

IF< IFGPL<  NjNTP  ).EQ.O  )  GO  TO  190 
DO  174  1=1,18 

TVP<  I  )=0.0 
DO  174  J=1 ,6 

174  TVP<  I  )=TVP(  I  )+BlA(  J ,  I  )*EPSDN<  J,N,NTP  ) 

K=3*II-2 
L=3#JJ-2 
M=3#KK-2 
DO  179  1=1,3 
J=I-i 

Pl<  K+J  )=P1<  K+J  >-TVP< I  ) 

Pl<  K+J+15  )=P1<  K+J+15  )-TVP<  1+9  ) 

Pl<  L+J  )=P1<  L+J  )-TVP< 1+3  ) 

Pl<  L+J+15  )=P1<  L+J+15 >-TVP< 1+12  ) 

Pl(  N+J  )=P1(  M+J  )-TUP<  1+6  ) 

179  Pl<  M+J+15  )=P1(  M+J+15  )-TVP< 1  +  15 ) 

190  CONTINUE 
RETURN 
END 

FUNCTION  NODEC I , J  ) 

COMMON/TD/ IMINl  100  >,1HAX< 100  ), JMINt 25  ), JMAX< 25  ), MAXI ,MAXJ»NMTL, NBC 
NQDE=0 

DO  100  JJ=1,J 
NSTAR'T=IMIN( JJ  ) 

NSTOP=IHAX( JJ  ) 

DO  100  II=NSTART ,NSTOP 
N0DE=N0DE+1 

IF<  JJ , EQ , J , AND .II.EQ.I  )  RETURN 
100  CONTINUE 
RETURN 
END 
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SUBROUTINE  POINTS 
INTEGER  CODE 

v  COMMON/BASIC/ ACELZ , ANGVEL , ANGACC ,  TREF ,  VOL , NUMNP  *  NUMEL,  NUMPC  , NUMSC , 

1NUMST 

COMMON /HATP/RO< 6  )»C< 12*16*6  )*EE( 16  )» AOrTS<  6 ) 

COMMON/NPDATA/R(  10  >,CQDE<10  ),XR<  10  ),Z<10  ),XZ<10  )» 

1NPNUM<  4,  8  ) ,  T  (  10  )>XT<  10  ) 

COMMON/ELDATA/BETA< 10  ),EPR<10  ),PR<4  ),SH<4  )»IX<8  ,5), 

1  IP(  4  )»JP<4  )» IS(  4  ).JS(4  )» ALPHA<  10  )»IT(4  )*JT<4  >» 

2ST( 4  ) 

COMMON/SOL VE/X(  883  ),Y<  888  ),TEM<  838  >»NUHTC  »  MBAND 

COMMON/TD/IMINC  100  ) » IMAX<  100), JMIN<  25 ),  JHAX<  25  ) » MAXI , MAX J , NMTL , NBC 
COMMON/PL ANE/NPP 

DIMENSION  AR<  4,  8>,AZ<  4,  8  )»MATRIL<  100*5  )»BLKANG<  100  )»BLKALF<  1 

100) 

DIMENSION  IBNG< 100  ),NBNG< 100  ) 

EQUIVALENCE  <  R<  1  ),AR  ),<  Z<  1  )» AZ  ) 

C  ESTABLISH  NODAL  POINT  INFORMATION 

C*  *********************************** 
NEL=0 
N0DSUM=O 
DO  100  J=1 ,MAXJ 
NSTART=IMIN( J  ) 

NSTOP=IMAX( J  ) 

DO  100  I=NSTART»NSTOP 
100  N0BSUM=N0DSUM+1 
NELSUM=0 
JJMAX=MAXJ-1 
DO  110  JJ=1 » JJMAX 
NSTOP=MINO(  IMAX(  JJ  )» IMAX<  JJ+1  )  )-l 
NSTART=MAXO< IMIN( JJ  ) , IMIN( JJ+1  ) ) 

DO  110  II=NSTART »NSTOP 
110  NELSUM=NELSUM+1 
NUMNP=N0DSUM 
NUMEL=NELSUM 
URITEC 15  )NUMEL » NUMNP 
DO  120  J=1 , MAXJ 
NSTART=IMIN( J  ) 

NSTOP=IMAX<  J  ) 

DO  120  I=NSTART  »NSTOP 
NPNUM<  I ,  J  )=NODE<  I » J  ) 

NP=NPNUM< I » J  ) 

R<  NP  )=AR<  If  J) 

120  Z<  NP  )=AZ<  I»J) 

C*  *********************************** 
C  READ  AND  ASSIGN  BOUNDARY  CONDITIONS 

C*  *********************************** 
C  INITIALIZE 

C*  *********************************** 
DO  130  1=1, NUMNP 
CODE<  I  )=0 

IF(  R<  I  ),EQ.O*  ♦  AND , NPP . EQ . 0  >  CODE(  I  )=1  ♦ 

XR<  I  )=0, 

XZ(  I  )=0 . 

XT(  I  )=0.0 
130  T(  I  >=0« 

IF<  NBC.EQ.O  )  GO  TO  210 
DO  200  IBC0N=1 ,NBC 

REAEK  5,1002)  II , 12, J1 , J2, ICN,RCON,ZCON,TCON 
DO  200  1=11,12 
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DO  200  J-Jl * J2 
NP=NPNUM< I,J  ) 
s  CODE(NP)=ICN 
XR(  NP  >=RCON 
XT( NP  )=TCON 
XZ(  NP  )=ZCON 
200  CONTINUE 
210  MPRINT=0 

WRITE<  15  X  CODE<  I  ),I  =  1,NUMNP> 

WRITE(  15  )(  XR<  I  )* 1=1 »NUMNP  ) 

WRITE<  15  )<  XT<  I  )» 1=1 »NUMNP  ) 

WRITE(  15  )(  XZ<  I  )» 1  =  1  *NUMNP  ) 

DO  230  J=1  ,MAXJ 
NSTART=IMIN( J  ) 

NSTOP=IHAX< J  > 

DO  230  I=NSTART »NSTQP 
NP=NPNUM< I»J  ) 

IF<  MPRINT .NE«0  )  GO  TO  220 
WR I TE<  6,2000) 

MPRINT=5? 

220  MPRINT=MPRINT-1 

230  WRITE<  6,2001  )  I, J,NP»CQDE(  NP  ),R<  NP  ),Z(  NP  ),XR(  NP  ),XZ(  NP  ),XT(  NP) 

C*  *********************************** 
C  ASSIGN  MATERIALS  IN  BLOCKS 

C*  *********************************** 
DO  300  M1=1,NUMEL 
300  IX<  Ml  ,5  )=0 

DO  310  IMTL=1 * NMTL 

READ  (5,1000)  MTL,(MATRIL<  IMTL,  IM  )» IM=2,5  ),BLKANG(  IMTL ),BLKALF(  IMT 
1L  ),  IBNG(  IMTL  ),NBNG(  IMTL  ) 

310  MATRIL( IMTL,1  )=MTL 

C*  *******************  **************** 
C  ESTABLISH  ELEMENT  INFORMATION 

C*  *********************************** 
J.JMAX=MAXJ-1 
N=0 
MTL=1 
KTL=1 

DO  440  JJ=1 , JJMAX 

NSTOP=MINO< IMAX< JJ  ), IMAX<  JJ+1  )  >-l 

NSTART=MAXO(  IMIN<  JJ )» IMIN< JJ+1  ) ) 

DO  440  II=NSTART,NSTOP 
NEL=NEL+1 

DO  400  IMTL=1 »NMTL 

IF< II *LT .MATRIL< IMTL, 2  ) )  GO  TO  400 

IF< II ♦GE.MATRIL< IMTL, 3  )  >  GO  TO  400 

IF<  JJ.LT »MATRIL< IMTL ,4 ) >  GO  TO  400 

IF( JJ.GE.MATRIL( IMTL, 5  ) )  GO  TO  400 

KAT=IMTL 

MAT=MATRIL(  IMTL,1 ) 

400  CONTINUE 

IF(KAT.EQ.KTL)  GO  TO  410 

KTL=KAT 

MTL=MAT 

GO  TO  420 

410  IF( II ,EQ,NSTART )  GO  TO  420 

IF<  JJ.NE, JJMAX, OR* II «NE« NSTOP  )  GO  TO  440 

M=NEL+1 

IANG=ICNG 

NANG-NCNG 
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GO  TO  421 

420  I--NPNUM<  IIiJJ) 

J=I  +  1 

K=NPNUM< II  +  1» JJtl  ) 

L=K-1 

M=NEL 

IX<  M»  1  )=I 

IX<  M» 2  )=J 

IX<  M>3  )=K 

IX<  M»4  )=L 

IX<  M»5  )=MTL 

BETA( M  )=BLKANG< KTL  ) 

ALPHA< M  >=BLKALF<  KTL  ) 

IANG=ICNG 
NANG=NCNG 
ICNG=IBNG<  KTL ) 

NCNG=NBNG( KTL ) 

421  NC=2 
430  N=N+1 

IF(M*LE«N)  GO  TO  440 
1X(  Nr  1  )=IX<  N-l » 1  )+l 
IX<  Nt 2  >=IX(  N-l  »2  )+l 
IX<  N*  3  )=IX(  N-l  »3  )+l 
IX(  N»4  )=IX( N-l  »4  HI 
IX(  Nf  5  )=IX<  N-l  »5  ) 

BETA(  N  )=BETA<  N-l  ) 

IF< IANG.EQ. 1  )  GO  TO  442 
ALPHA<  N  )=ALPHA(  N-l  ) 

GO  TO  443 

442  CONTINUE 
IF(NC.GT.NANG)  GO  TO  444 
ALPHA(  N  )=ALPHA( N-l  ) 

GO  TO  443 
444  NC=1 

ALPHA(  N  )=-ALPHA< N-l  ) 

443  CONTINUE 
NC=NC+1 

IF(M.GT.N)  GO  TO  430 
440  CONTINUE 

IF< NUMNP.GT.2000  )  URITE< 6 » 2002 > 

C*  *********************************** 
C  SET  NODAL  POINT  TEMPERATURE  TO  REFERENCE  TEMPERATURE 
C*  **************  ********************* 
IF<  NUMTC.NE.O )  RETURN 
DO  500  N=1 *NUMNP 
500  T<  N  )=TREF 

1000  FORMAT  < 5I5»2F10.0»2I5 ) 

1002  FORMAT<  4I5»I10»3F10»0) 

2000  FORMAT  < 128H1  I  J  NP  TYPE  R-ORDINATE  Z-ORDINA 

1TE  R  LOAD  OR  DISPLACEMENT  Z  LOAD  OR  DISPLACEMENT  T  LOAD  OR  DISP 
2LACEMENT  ) 

2001  FORMAT  < 2I5> I4» I12»F13.6 f F14 . 6»E26 .7,E24.7>E24.7  ) 

2002  FORMAT  < 35H  BAD  INPUT  -  TOO  MANY  NODAL  POINTS) 

RETURN 

END 

SUBROUTINE  QUAD 
INTEGER  CODE 

REAL  NUSN >  NUTN >  NUTS i NUNS » NUNT » NUST 
DIMENSION  DUMMY<6,A)»DUMMY1<6>6) 

C0MM0N/PLAS/ALFA<6*  4»8  )»SIGYLD( 7»6»3  )f IFGPL(  4,8) 
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COMMON /ARG1/SIG1<  18>,EPS1<  18)»BEP3P;  12)»CEP£p: 
COMMON/NXDAT  A/NTP  >  NTS  *  NTGT3 » GTS1G?  24  >  24  *  S  ) 
COMMON /BAS I C/ACELZ , ANGVEL , ANGACC » TREE ? VOL »  NUMNP 


6  ■■  6  ) 

, NUiTEL  > 


NLiilPC » KUM5C » 


1NUMST 


COMMON/NXC1UAD/AR1 

COMMON /NON AX I /S 1< 30 >30  )>P1< 30  )» THETA? BS1< 6.30  ) 


COMMON/MATP/RO(  6  )?E<  12»16,6)>EE(  16  ) >  AOFTS(  6  )  • 


COMMON/NPDATA/R<  10  >, CODE?  10  >,XR<  10  )>Z<10  ),XZ<10  ), 

1NPNUM(  4 >  8  )>  T<  10  )>XT(  10  ) 

COMMON/ELDATA/BETA<  10  )>EPR(10  )>PR<  4  )#SH<  4  ),IX<8  >5), 

1  IP(  4  )>  JP(  4  ),  IS<  4  )>  JS(  4  )>  ALPHA?  1  0  )»IT<4  )»JT<4  )> 

2ST< 4  ) 

COMMON/ ARG/RRR(  5  ),ZZZ<  5  ) »  RR<  4  >,ZZ<  4  ),S<  15.15  )>P<  15  )>TT<  6  )> 

1H(  6>  15  )*CRZ(  6>6  )»XI<  10  ) >  ANGLE?  4  )>SIG(  18  )>EPS(  13)>N 
COMMON/INCR/NOLINC»  NOL  > INERT  >  NUMMAT  >  SIGTQT<  12.  4.8) 


COMMON/RESULT/BS(  6>15  ),D<  6 >6  ),C(  6>6  ),AR,BE«  6>9  )>CNS(  6,6  ) 
COMMON/PLANE/NPP 

C0MM0N/DUM1/S1TEM( 3 >30  ),S1T< 24>24  )>TS(  6>24) 

DIMENSION  S2T<  24  >  6  ) 

DIMENSION  BS1T<  6>3  )  ,P1T<3>  > PITT?  24) 

1 1  =  IX(  N>  1  ) 

J1  =  IX(N>2) 

K1=IX<  N>3  > 

L1=IX<N,4) 

MTYPE=IX(N»5) 

IX(  N>5  )=-IX(  N>5  ) 


c* 

* 

*********** 

****** 

* 

*  * 

* 

* 

*********:( 

k  ^  JjC 

c 

c* 

* 

INTERPOLATE  MATERIAL 
*********** 

PROPERTIES 

****** 

* 

*  * 

* 

* 

*********  3 

k  *  * 

DO  100  1=1,12 
100  EE<  I  )=E<  1,I+1,MTYPE  > 

DO  110  1  =  1 >  6 
DO  110  .J=l  ,6 
CNS(  I  >  J  )=0.00 
C<  I, J  )=0.00 
110  D<  I  >  J  )=0.00 

C*  *******************  *  4  *  *  *  *  *********  * 
C  FORM  STRESS-STRAIN  RELATIONSHIP  IN  N-S-T  SYSTEM 

C*  ********************  *  *****  *  *******  * 
NUNS=EE( 4  ) 

NUNT=EE( 5  ) 

NUST=EE< 6  ) 

NUSN=<  EE< 2  )*NUNS  )/EE( 1  ) 

NUTN=<  EE(  3  )*NUNT  )/EE(  1  ) 

NUTS=(  EE<  3  >*NUST  )/EE(  2  ) 

•DI  V=1 . 00-NUNS*NUSN-NUST*NUTS-NUNT*NUTN-NUSN*NUNT*NUTS 
1-NUNS*NUTN*NUST 

CNS<  1>1  >=EEC1  >*<  1.00-NUST*NUTS)/DIV 
CNS(  1 ,2  >=EE<  2  )*<  NUNS+NUNT*NUTS  )/DIV 
CNS( 1 >3  )=EE( 3 )*<  NUNT+NUNS*NUST )/DIV 
CNS<  2,1  )=CNS(  1  >  2  ) 

CNS(  2,2  )=EE(  2  )*<  1 .OQ~NUNT*NUTN  )/DIV 
CNS( 2,3  )=EE( 3  )*<  NU3T+NUSN*NUNT  )/DIV 
CNS<  3,1  )=CNS<  1,3  ) 

CNS(  3,2  )=CNS<  2,3  ) 

CNS( 3,3  )=EE( 3  )*< 1 .00-NUNS*NUSN )/DIV 
CNS<  4,4  )=EE(  7  ) 

CNS(  5,5  )=EE(  8  ) 

CNS<  6,6  )=EE<  9  ) 

DO  162  1=1,6 
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DO  162  J=l»o 
162  CEPSPt  I»J>=0.0 

IF  (  IFGF'L<  N.NTP  ).NE.O  )CALL  ELPLSSC  MTYPE  > 

SET  UP  STRAIN  TRANSFORM  TO  N-S-T  SYSTEM 
SINA=SIN<  ALPHAC N  )  ) 

COSA=CQS( ALPHA(  N  ) ) 

S2=SINA**2 
C2-C0SA:K^2 
SC=SINA*CCSA 
D<  1  *  1  )=C2 
D<  If  3  >=S2 
D<  1  f  6  )=-SC 
D<  2f  1  )=S2 
EK  2f 3  )=C2 
D<  2f6  )=SC 
EK  3f 2  )=1 . 00  . 

D<  4f 1  )=2.00*SC 
EK  4f  3  >=-2.00*SC 
D< 4f6  )=C2-S2 
D<  5f  4  )=SINA 
EK  5f5  )=CQSA 
EK  6f4  )=CQSA 
EK  6f  5  )=-SINA 

C  SET  UP  STRAIN  TRANSFORMATION  TO  R-Z-T  SYSTEM 
SINB=SIN<  BETA(N)) 

COSB=COS< BETA( N  )  ) 

S2=SINB**2 
C2=C0SB**2 
SC=SINB*CQSB 
C(  If  1  )=S2 
C<  1  f  2  )=C2 
C<  lf4)=SC 
C(  2  f  1  )=C2 
C(  2f  2  )=S2 
C(  2f  4  )=-SC 
C(  3f  3  )=1 .00 
C(4f  1  )=-2,00*SC 
C<4f2>=2.00*SC 
C<  4f  4  )=S2-C2 
C<  5f  5  )=SINB 
C(  5f 6  )=-COSB 
C(  6f  5  )=COSB 
C(  6f 6  )=SINB 

IF  < IFGPL<  N f  NTP  ) . NE . 0  )CALL  ELPLSS( MTYPE  ) 

C  CALCULATE  CRZ  MATRIX 
DO  120  1=1 f 6 
DO  120  J=1 f 6 
DUMMY(  If  J  >=0.00 
DO  120  K=1 f 6 

1 20  DUMMY(  I  f  J  >=DUMM  Y<  I  f  J  )+D<  I  f  K  )*C(  K  f  J  > 

DO  130  1  =  1  f  6 
DO  130  J=1 f 6 
DUMMY  1<  If  J  >=0.00 
DO  130  K=1 f 6 

130  DUMMYK  I  f  J  >=DUMMY1<  I  f  J  )+CNS(  I  f  K  >*DUMMY<  Kf J  > 
DO  140  1  =  1  f  6 
DO  140  J=1 f 6 
DUMMY<  I  f  J  >=0 . 00 
DO  140  N=1  f  6 

1 40  DUMMY(  I .  J  )=DUMMY(  I  f  J  >+D<  K » I  >*DUMMY  1(  K » J  > 
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HO  160  1  =  1  f  6 

DO  150  J=1  f  6 
CRZ<  I»J  >=0,00 
DO  150  K=1 » 6 

150  CRZC  If  J  )=CRZ(  If  J  )+C<  K » I  )$BUHMY(  K»  J  ) 

'  TT<  I  )=0«00 

no  160  M=l»6 

P<  M  >=0.00 

no  161  1 1  =  1 » 3 

IF<  AOFTS(  MTYPE  ).EQ.  1 .  )  P<  M  )=CNSC  H»  II  >*EEC  III?) 
161  P<  H  )=P<  M  >+(  T<  N  )-TREF  >*CNSC  M  f  1 1  >*EEC  1 1  f?  > 

00  160  K=1  »6 

160  TT<  I  >=TT(  I  >+CCKf  I  >*DC  MtK  )*PC  M  > 

C 

C  FORM  QUADRILATERAL  STIFFNESS  MATRIX 

RRR(  5  )=(  RC  II  >+R(  J1  )+R(Kl  HRCLl  >>/4. 

ZZZ<  5  >=<  ZCI1  >+Z(  J1  >+Z<  K1  >+Z<  LI  )  )/4 . 
no  200  M=1 f 4 
MM=IX(  Nf M ) 

IF<  NPP.NE.O  )  GO  TO  190 

IF< R<  MM  >.EQ.O. .AND. CODEC  MM  > .EQ . 0 . >CODE(  MM >=1 . 
190  RRRC  M  >=R<  MM  ) 

200  ZZZC  M  )=Z<  MM  ) 
no  220  11=1 t 15 
Pl(  II  >=0.0 
PIC  11+15)  =0.0 
PC  II  )=0.00 
DO  220  JJ=1 » 15 
220  SC  IIfJJ>=0.00 
V0L=0. 

DO  90  1=1 t 6 
DO  90  J=1 » 15 
BS1C  If  J  )=0.0 
BSK  If  J+15  )  =  0.0 

90  BSC  I  f  J  )=0.00 
AR=0 . 00 

240  CALL  TRISTFC4flf5> 

CALL  TRISTFC If2f5) 

CALL  TRISTFC 2f3f 5) 

CALL  TRISTFC 3f4f 5 ) 

DO  91  1=1 f  6 
DO  91  J=lf 15 

91  BSC  IfJ)=BSC  If  J)/AR 

DO  300  1=1 f 30 
DO  300  J=1 f 30 
300  SIC  If  J  >=0.0 
AR1  =0.0 

CALL  NAXSTFC4flf5> 

CALL  NAXSTFC If2f5) 

CALL  NAXSTFC 2f3f 5) 

CALL  NAXSTFC 3f4f 5) 

DO  310  1=1 f 6 
DO  310  J=1 f 30 

310  BS1C  If  J  >=  BS1C  If  J  >/ARl 
DO  320  1=1 f 6 
DO  320  J=1 f 3 

320  BSITCIfJ)  =  BSlCIfJ+12> 

DO  325  1  =  1  f  6 
DO  325  J=1 » 12 

325  BS1C  I  f  J+12  >  =  BSlCIfJ  +  15> 
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DO  330  1=1 *6 
DO  330  J=1 » 3 

330  BS1(  I » J+24  )  =  BSlT(ItJ) 

DO  340  1=1*3 

340  PIT  (  I  )  =  Pl<  1  +  12) 

DO  341  1=1*12 

341  Pl< 1  +  12  >  =  Pl< 1  +  15  ) 

DO  342  1=1*3 

342  Pl< 1+24  )  =  P1T< I  ) 

DO  149  1=1*3 

DO  149  J=1 *30 

149  S1TEM(  I » J  )  =  SHI  +  12.J) 

DO  151  1=1*12 
DO  151  J=i *30 

151  Sl< 1+12* J  )  =  Sl( 1+15* J  ) 

DO  152  1=1*3 

DO  152  J=l*30 

152  Sl(  1+24* J  )  =  S1TEM(I,J) 

DO  153  1=1,3 

DO  153  J=1 *30 

153  S1TEM< I  * J  )  =  Sl< J*I  +  12  ) 

DO  154  J=1 *  12 

DO  154  1=1,30 

154  Sl( I » J+12  )  =  Sl( I  * J+15  ) 

DO  155  1=1,3 

DO  155  J=1 ,30 

155  Sl<  J *1+24  )  =  S1TEM( I » J  ) 

DO  251  1=1,6 
DO  251  J=1 *24 

251  TS< I , J  )  =  0.0 

DO  252  1=1,3 
DO  252  J=1  *4 
TS< I,I+< J-l >*3  )  =  0.250 

252  TS<  I+3,I+12+<  J-l  >*3)  =  0.250 

DO  253  1=1,24 
DO  253  J=1 *24 
S1T< I* J  )  =  0.00 
DO  253  K=1 *6 

253  S1T<  I*  J  )  =  S1T<  Z  *  J  )  +  S1<I»24+K)*TS(K,J) 

DO  254  1=1,24 

DO  254  J=1 *24 

254  Sl<  I  *  J  )  =  S1(I»J)  +  S1T<I,J)  +  S1T<  J*  I  > 

DO  255  1=1*24 

DO  255  J=1 ,6 
S2T(  I* J  )  =  0.0 
DO  255  K=1 *6 

255  S2T( I  * J  )  =  S2T< I* J  )  +TS< K, I  )*S1< K+24* J+24 > 

DO  256  1=1*24 

DO  256  J=1 *24 

51T( I* J  )  =  0.0 

DO  256  K=1 ,6 

256.  S1T<  I » J  )  =  S1T(  Z* J  )  +  S2T(I»K)*TS<K*J> 

DO  257  1=1,24 

DO  257  J  =1,24 

257  S1(I*J)  =S1<  I » J  )+SlT<  I » J  ) 

DO  258  1=1*24 

P1TT<  I  )=0.0 
DO  258  K*1 *6 

258  P1TT<  I  )=P1TT<  I  >+TS<  K»  I  >*PK  K+24  ) 

DO  259  1=1,24 
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no  non  non  non 


259  Pl<  I  )=P1<  I  ) fPITT (  I  ) 


RETURN 

END 

SUBROUTINE  SOLV 

COMMON/ELDATA/BETA< 10  >»EPR<  10  >»PR(4  >»SH<4  MX(  3  »5>» 

1IP<  4  )»JP<  4  ) » I  S(  4  )» JS<  4  )  i  ALPHA;  10  )»IT(4  )*JT<4  )» 


2ST<  4  ) 

COMMON/BASIC/ACELZ , ANGVEL *  ANGACC *  TREE * VOL * NUKNP 
1NUMST 

C0MM0N/N0NAXI/SK  30*30  >»P1<  30  >* THETA* BS1<  6*30  ) 
COMMON/NXD AT A/NTP *  NTS  *  NTOTS *  GTS 1G(  24  *  24  *  S ) 
C0MM0N/S0LVE/B<  72)>A<  72*36)*  NUMBLK  *  ft BAND 


NUMEL* NUNPC  *  HUM SC* 


MM=MBAND 

NN=36 

NL=NN+1 

NH=NN+NN 

REWIND  1 

REWIND  2 

NB=0 

GO  TO  150 


C* 

* 

*  *  *  * 

tj.  ^  ^  ^ 

*  *  *  *  * 

* 

* 

* 

* 

* 

* 

* 

*  * 

* 

* 

*  i 

t  *  *  * 

* 

*  *  * 

c 

c* 

* 

REDUCE 
*  *  #  * 

EQUATIONS 
#  *  #  *  # 

BY  BLOCKS 
#  #  *  *  * 

* 

* 

* 

* 

* 

* 

* 

*  * 

* 

4 

t  i 

h  *  % 

* 

*  *  * 

1.  SHIFT  BLOCK  OF  EQUATIONS 


100  NB=NB+1 

DO  125  N=1 » NN 
NM=NN+N 
B<  N  )=B<  NM  ) 

B(NM  >=0.00 
DO  125  M=1 »MM 
A(N,M)=A<  NM»M> 
125  A<  NM»M  )=0.00 


2.  READ  NEXT  BLOCK  OF  EQUATIONS  INTO  CORE 


IF<  NUMBLK.EQ.NB  )  GO  TO  200 
150  REAIK2)  (B(N>»<A(N*M)fM=l*MM)»  N=NL  *  NH  ) 
IF(NB.EQ.O)  GO  TO  100 

3.  REDUCE  BLOCK  OF  EQUATIONS 


200  DO  300  N=1»NN 

IF<  A( N» 1  KEQ.0.00  )  GO  TO  300 
B(  N  )=B<  N  )/A<  N»  1  ) 

DO  275  L=2»MM 

IF(A(N*L).EQ.0.00)  GO  TO  275 
C=A(N»L)/A(N»1  ) 

I=N+L-1 

J=0 

DO  250  K=L»MM 
•J=JF1 

250  A(  I*.J)=A<  I*J)-C*A<N,K) 

B<  I  )=B(  I  )-A<  N»L  >*B<  N) 

A<  N»L  >=C 
275  CONTINUE 
300  CONTINUE 

4.  WRITE  BLOCK  OF  REDUCED  EQUATIONS  ON  FORTRAN  UNIT  1 
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c 

I F<  NUMBLK .  EQ . NB )  GO  TO  400 

WRITE  (1)  <B<N)»(A<N>M)*M=2»MM)*N=1»NN) 

GO  TO  100 

C*  *********************************** 
C  BACK-SUBSTITUTION 

C*  *********************************** 
400  DO  450  M=1 f NN 
N=NN+1-M 
DO  425  K=2»MM 
L=N+K-1 

425  B<  N  )=B<  N  )-A(  N»K  )*B<  L  ) 

NM=N+NN 
B<  NH  >=B<  N  ) 

450  A<  NM»NB  )=B(  N  > 

NB=NB-1 

IF<  NB *EQ *0  )  GO  TO  500 
BACKSPACE  1 

READ  (1)  (  B(  N  )»<  A<  N»M  )»M=2»MM  )»N=1»NN  ) 

BACKSPACE  1 
GO  TO  400 

C*  *********************************** 
C  ORDER  FORMER  UNKNOWNS  IN  B  ARRAY 

C*  *********************************** 

500  K=0 

DO  600  NB=1 »NUMBLK 
DO  600  N=1»NN 
NM=N+NN 
K=K+1 

600  B(  K  )=A<  NM»  NB  ) 

C*  *********************************** 
C  WRITE  SOLUTION 

C*  *********************************** 
NN12  =  3*NUMNP 
1500  FORMAT(  "  "  »5I10) 

WRITE<  26  )  <  B(  I  )»I  =  1»NN12) 

WRITEC 26 X < IX< I » J ) » J=1 ,5 >» 1=1 .NUMEL ) 

MPRINT=0 

DO  710  N=1 »NUMNP 

IF<  MPRINT . NE»0  )  GO  TO  700 

WRITE  <  6 » 2000  ) 

MPRINT=5? 

700  MPRINT=MPRINT-1 

710  WRITE  ( 6 *2001  )  N,B< 3*N-2  )»B( 3*N-1  )»B(  3*N  ) 

2000  FORMAT  ( 13H1  NODAL  P0INT>18X»2HUR»18Xr2HUZ»18X»2HUT ) 

2001  FORMAT  (I13»3E20.7) 

RETURN 

END 

SUBROUTINE  STIFF 
INTEGER  CODE 

COMMON/RESULT /BS<  6fl5)»D(6»6)tC(6»6  )» AR»BB<  6t  9  )» CNS<  6.6) 

COMMON/ ARGl/SIGK  18)»EPS1<  18)rBEPSP<  12 )» CEPSPt 6*6  ) 
COMMON/BASIC/ACELZ  *  ANGVEL  *  ANGACC » TREF  *  VOL  »NUMNP » NUMEL  »NUM  PC » NUNSC  * 
1NUMST 

COMMON/ELDATA/BETA<  10  ),EPR(  10  >»PR<4  >*SH<4  >»IX<8  *5)» 

1IP(4  )*JP(4  )*  IS<  4  )» JS<  4  )»ALPHA<  10  )»IT<  4  )*JT(4  )* 

2ST<  4  ) 

COMMON/NPDATA/R<  10  )*CODE<  10  >.XR<  10  )*Z(10  )»XZ<10  )* 

1NPNUM<  4*  8  )»T(  10  )*XT<10  ) 

COMMON/SOL VE/B<  72).A<  72.36  )>NUMBLK.MBAND 
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COMMON/NXDAT A/NTP , NTS , NTOTS , GTS1G<  24,24,8) 

C0MMQN/ANS4/FT<  24,4  ),GTS1U<  24  ),GTS1UT<  24,4) 

CQMMON/ARG/RRR(  5  ) , ZZZ<  5  ) , RR(  4  ) , ZZ(  4  ) » S(  15,15), P<  15),TT(6), 

1H(  6, 15  ),CRZ(  6,6  )»XI<  10  ),  ANGLE<  4  )»SIG<  18)»EPS<  18),  N 
-  COMMON/NONAXI/SK  30*30  ),P1<  30), THETA, BS1<  6,30) 

COMMON/PLANE/NPP 

C0MM0N/ANS2/GTP1<  24 ),G<  24,24 ),GTS1<  24,24  )»GTS1GE<  24,24  ) 
C0MM0N/BUM1/S1TEM(  3,30  ),S1T<  24,24  ),TS(  6,24  ) 

COMMON/RATE/DKPR , SIGPR , BVR , EUR , PSRATE(  10,8 ),  NRATE 
DIMENSION  LM<  4  ),S2(  12,3  ),S3(  3»12),S4<  3,3),S5<  12,3  ),S6(  12,12  ) 

C*  ************************  *********** 
C  INITIALIZATION 

NRATE=1 
REWIND  2 
REWIND  3 
NB=12 
ND=3*NB 
ND2=2*ND 
ST0P=0. 

NUMBLK=0 
DO  100  N=1 ,ND2 
B<N>=0.00 
DO  100  M=1»ND 
100  A(  N,M  )=0 .00 
DO  50  1=1,24 
FT<I,NTP)  =  0.0 
GTS1UT<  I ,NTP  )=0.0 
DO  50  J=1 ,24 

50  GTS1G< I , J,NTP  )  =  0.0 

Ct  *********************************** 

C  FORM  STIFFNESS  MATRIX  IN  BLOCKS 

Ct  *********************************** 
200  NUMBLK=NUMBLK+1 
NH=NB*( NUMBLK+1  ) 

NM=NH-NB 

NL=NM-NB+1 

KSHIFT=3*NL-3 

DO  340  N=1 ,NUMEL 

IF<  IX<  N,5).LE.0)  GO  TO  340 

DO  210  1=1,4 

IF<  IX<  N,I).LT.NL)  GO  TO  210 
IF<  IX(N,I  ).LE.NM)  GO  TO  220 
210  CONTINUE 
GO  TO  340 
220  CALL  QUAD 

IF< VOL.GT.O.  )  GO  TO  230 
WRITEC6»2000 )  N 
STQP=1 » 

230  IF<  IX(N,3  ).EQ.IX<  N,4  ) )  GO  TO  300 
DO  231  11=1,3 

DO  231  JJ=1»3 

231  S4<II,JJ)=S<  11+12,  JJ+12) 

CALL  SYMINV<  S4,3 ) 

DO  232  11=1,12 
DO  232  JJ=1 ,3 

232  S2<  II,  JJ  >=S<  II,  JJ+12) 

DO  233  11=1,3 

DO  233  JJ«1 , 12 

233  S3<  II, JJ)=S<  11+12, JJ) 

DO  240  1=1,12 
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no  240  J=1 ,3 
S5<  I.J  )r=0*00 
no  240  K=1 , 3 

240  S5<  I » J  >  =  S5(  I » J  )  +  S2(  I » K  )  *  S4<K,J> 
no  241  1=1,12 

no  241  J=1 , 12 
S6<  I»  J  >=0.00 
DO  241  N=i ,3 

241  S6<I,J>  =  S6<  I » J  )  +  S5(  I »K )  *  S3<K,J) 

DO  234  I 1=1 » 12 

DO  234  JJ=1 *3 

234  P<  II  )=P<  II  )-S5<  II,  JJ  >*P<  JJ+12) 

DO  235  11=1.12 

DO  235  JJ=1 » 12 

235  S<  II.JJ  >=S<  II,JJ  >-S6<  II,  JJ  > 

DO  259  1=1.24 

DO  259  J=1 ,24 

259  G< I.J )  =  0.0 
DO  260  K= 1 » 4 
DO  260  1=1,3 

G< K*3-3+I, 1*4-3 )  =1.0 
G<K*3-3+1, 1*4-2)  =  RRR<  K  ) 

G< K*3-3+I, 1*4-1 >  =  ZZZ( K  ) 

260  G(K*3-3+I,I*4  )  =  ZZZ<  K  )  *RRR<  K  > 

DO  262  1=1,12 
I»0  262  J=1 ,12 
262  G<  1+12, J+12  )  =  G( I , J  ) 

NTP20  =  21 

WRITE<  NTP20  )  (<  CRZ< I, J >,J=1,6 >.1=1,6 > 

WRITE( NTP20 X  <  BS1< I, J  ), J=l»30  >,1=1,6 ) 

WRITE<NTP20X<G<  I.J  ),J=1,24  >,1  =  1,24) 

WRITE< NTP20 X  <  CEPSPC I ,  J  >,  J=1 ,6  ),  1=1 , 6  ) 

WRITE<  NTP20  X  <  CNS<  I,J),J=1,6),I=1,6> 

WRITE< NTP20  X  <  D< I , J  ) , J= 1 , 6  ) , 1=1 , 6 ) 

WRI TE<  NTP20  X  <  C< I , J  ) , J= 1 , 6  ) , 1= 1 , 6  > 

DO  280  1=1,24 
GTP1<  I  >=0.0 
DO  280  K= 1,24 
GTS1< I »K  )  =  0.0 

GTP1<  I  )=  GTPKIH  G<K,I)*P1<K> 

DO  280  J=1 ,24 

280  GTSKI.K)  =  GTSKI.K)  +  G(J,I)  *  S1(J,K) 

WRITE<  3)  <<GTS1<  I.J  >,J=  1,24  >,1  =  1,24) 
no  281  1=1,24 

FT<  I  ,NTP  >  =FT<  I  ,NTP  )  +  GTP1< I > 

DO  281  J=1 ,24 
GTS1GE< I , J >  =  0.0 
DO  281  K=1 ,24 

281  GTS1GE<  I ,  J  )  =  GTS1GE<  I,  J )+  GTSKI.K)  *G<K,J> 

DO  282  1=1,24 
DO  282  J=1 ,24 

282  GTS1G< I » J»NTP )  =  GTS1G< I , J.NTP >  +  GTS1GE(  I, J  ) 

C*  *********************************&$ 
C  ADD  ELEMENT  STIFFNESS  MATRIX  TO  BODY  STIFFNESS  MATRIX 
C*  *********************************** 
300  DO  310  1=1,4 
310  LM<  I  >=3*IX<  N,I  )-3 
DO  330  1=1,4 
DO  330  K*1 ,3 
I  I=LM<  I  HK-KSHIFT 


KK=3*I-3+K 

B<  II  )=B<  II  KP<  KK  ) 

DO  330  3=1  >4 

DO  330  L=l>3 

33=LM<  J  WL--I I  +  l-KSHIFT 

LL=3fc3~3TL 

IF(  J.J.LE.O  )  GO  TO  330 
IF(ND.GE.JJ)  GO  TO  320 
WRITE( 6 t 2001  )  N 
STOF-1. 

GO  TO  340 

320  A<  1 1 ,  J.J  )=A(  IlfJJ  )+S<  KK»LL  ) 

330  CONTINUE 
340  CONTINUE 

C*' ft********************************** 

C  ADD  CONCENTRATED  FORCES 

C*  ******************#*: %.%%%%**%*%*<*%%* 
DO  400  N=NL»NM 
IF< N.GT.NUHNP  )  GO  TO  500 
K=3#N-KSHIFT 
B(  K  )=B<  K  )+XT (  N  ) 

B<  K-l  >=B<  K-l  )+XZ<  N  ) 

400  B<  K-2  )=B(  K-2  )+XR<  N  ) 

C  ******4*********************:*#****** 

C  ADD  PRESSURE  BOUNDARY  CONDITIONS 

500  IF<  NUMPC.EQ.O )  GO  TO  600 
DO  540  L=1 »NUMPC 
I  =  IP(  L  ) 

J=3P<  L  ) 

PP=PR(  L  )/6  * 

BR=(  R<  J  )-R(  I  )  )*PP 
DZ=<  Z<  I  )-Z<  J  )  >*PP 
RX=2.*R<  I  )+R<  J  ) 

ZX=R<  I  )+2.*R(  J) 

II=3*I-KSHIFT-1 

33=3*J-KSHIFT-1 

IF< II .LE.O.OR. II .GT.ND )  GO  TO  520 
SINA=0  • 

C0SA=1. 

510  B< II-l  )=B< II-l  )+RX*<  COSA*DZ+SINA*DR ) 

GR=RX*<  COSA*DZ+SINA*DR )*THETA/2.0 
FT(  1 » NTP  )  =  FT(  1  *NTP  )+GR 
FT(  2»NTP  )  =  FTC  2, NTP)  +R<  I  )*GR 
FT(  3»NTP  )  =  FT(  3*NTP  )  +ZC  I  )*GR 
FT<4»NTP)  =  FT<  4»NTP  )  +  R(  I  >*Z(  I  )*GR 
FTC  14»NTP  )=  FT(  14  f  NTP  >  +R<  I  )*GR 
FT< 13»NTP  )  =  FTC 13 » NTP  )  +GR 

FT< 15»NTP  >  =  FT< 15»NTP  )+Z( I  )*GR 
FT(  16 » NTP )  =  FT< 16»NTP )  +Z(  I  >*R<  I  )#GR 
B(  II )=B< II  )-RX*<  SINA*DZ-COSA*DR  ) 

GZ=-RX*<  SINA*DZ-COSA*DR )  *THETA/2 . 0 
FT<  5»NTP  )=FT<  5*NTP  )+GZ 
FT(  6rNTP  )=FT<  6>NTP  )+R(  I  )*GZ 
FT(  7 »NTP  )=FT<  7 »NTP  )+  Z< I >*GZ 
FT< 8»NTP  )=FT( 8»NTP  )+Z< I  )*R< I  )*GZ 
FT<  17, NTP  )=FT<  17, NTP  HGZ 
FTC  18, NTP  )=FT(  18, NTP  HR<  I  >*GZ 
FTC  19,NTP)=FT<  1?,NTP)+Z<  I  >*GZ 
FTC  20, NTP  )=FTC  20,NTPHZC  I  )*R<  I  )*GZ 


-  V;,*' 
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520  IF<  JJ  *LE ♦ 0 »0R  *JJ  «GT * ND )  GO  TO  540 
SINA=0» 

CQSA=1. 

530  E'.<  JJ-l  )=B<  JJ-l  >+ZX*<  COSA*BZ-fSINA*DR  ) 

GR=  ZX  *C C0SA*DZ+S1NA*BR)  *THETA/2.0 
FTC  1 .  NTP  >=FTC  1  »NTP  )+GR 
FTC 2. NTP >=FTC 2.NTP  HRC J  )*GR 
FTC3.NTP)=FTC3.NTPHZC  J  >*GR 
FT<  4.  NTP  >=FTC  4»NTP  )+Z(  J  >*RC  J  )*GR 
FT<  13 .NTP  )-FTC  13. NTP  )+GR 
FTC  14. NTP  )=FTC  14. NTP  )+RC  J  )*GR 
FTC  15. NTP  )=FTC  15. NTP  HZC  J  )*GR 
FTC  16.NTP  )=FTC  16.NTP  )+2<  J  >*RC  J  )*GR 
BC  JJ  >=BC  JJ  >-ZX*C  SINA*DZ-COSA*DR  ) 

GZ=  -ZX*C SINA*BZ-COSA*DR>  *THETA/2.0 
FTC  5.NTP  )=FTC  5.NTP  >+GZ 
FTC  6.  NTP  >=FTC  6.NTP  )+RC  J  >*GZ 
FTC  7 , NTP  )=FTC  7 » NTP  )  +ZC  J  >*GZ 

FTC  8 .  NTP  )=FTC  8 » NTP  )+ZC  J )*RC  J  >*GZ  * 

FTC  17. NTP  )=FTC  17.NTPHGZ 
FTC  18 . NTP  )=FTC  18.NTP  HRC  J  )*GZ 
FTC  1?.NTP)=FTC  19.NTPHZC  J  >*GZ 
FTC  20. NTP  )=FT( 20.NTP  HZC  J  HRC  J  >*GZ 
540  CONTINUE 

1100  FORMATC "  " .12E10.3) 

C*  *********************************** 

C  ADD  SHEAR  BOUNDARY  CONDITIONS 

C*  *********************************** 
600  IFC  NUMSC . Ed . 0  )  GO  TO  701 
DO  640  L=l» NUMSC 
I=ISCL) 

J=JSC  L> 

5S=SHC  L  )/6 . 

DZ=C  ZC  I  )-ZC  J  )  )*SS 
DR=<  R(  J  )-R<  I  )  >*SS 
RX=2.*RC  I  )+RC  J  ) 

ZX=RC  I  H2.*RC  J  ) 

II=3*I-KSHIFT-1 

JJ=3*J-KSHIFT-1 

IFC II.LE.O.OR.II.GT.ND)  GO  TO  620 
SINA=0» 

C0SA=1 • 

610  BC  I.T-1  )=B(  II-l  HRX*C  SINA*DZ+COSA*DR  ) 

GR=  RX*C SINA*DZ+COSA*DR)  *THETA/2.0 
FTC  If  NTP  )=FTC  l.NTPHGR 
FTC  2.  NTP  )=FTC  2.NTP  HRC  I  >*GR 
FTC  3,  NTP)  =FTC3»NTPHZ<  I  >*GR 
FTC  4» NTP  )=FTC4»NTP)+ZC  I  )*RC  I  >*GR 
FTC  13.NTP  )=FTC  13.NTPHGR 
FTC  14. NTP  )=FTC  14. NTP  >+RC  I  )*GR 
FTC  15.NTP  )=FTC  15.NTP  )+ZC  I  )*GR 
FTC 16.NTP )=FTC 16. NTP HZC I  HRC I  HGR 
BC II  )=BC II  >-RX*<  COSA*DZ-SINA*DR  ) 

GZ=  -RX*CCOSA*DZ-SINA*DR>  *THETA/2.0 
FTC  5f NTP  >=FTC  5.NTP )+GZ 
FTC  6.NTP  )=FTC  6, NTP  >+RC  I  HGZ 
FTC  7  .NTP  )=FTC  7  »NTP )+ZC I HGZ 
FTC  8» NTP  )=FTC  8.  NTP  HZC  I  HRC  I  HGZ 
FTC 17.NTP  )=FTC 17. NTP HGZ 
FTC  18.NTP  )=FTC  18. NTP  HRC  I  HGZ 
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FT( 19»NTP  )=FTC 19, NTP  )+Z(  I  )*GZ 
FT< 20, NTP  )=FT(  20, NTP  KZC I  )*RC I >*GZ 
620  IF<  JJ*LE.0*0R. JJ .GT»ND )  GO  TO  640 
SINA=0. 

C0SA=1. 

630  Ef<  JJ-1  >=B<  JJ-1  HZX*C  SINA*HZ+COSA*DR  ) 

GR=  ZX#(  SINA$DZ+COSA#BR  )$THETA/2 . 0 
FT< 1  »NTP  )=FT< 1  ,NTP  )  fGR 
FT<  2»NTP  )=FT<  2 ,  NTP  >+RC  J  )*GR 
FT(  3»NTP  )=FTC  3,  NTP  HZC  J  )*GR 
FT<  4,  NTP  )=FT<  4, NTP  HZC  J  )*R<  J  )*GR 
FT<  13»NTP  )=FT<  13»NTP  HGR 
FT<  14»NTP  )=FT(  14, NTP  HRC  J  HGR 
FT<  15»NTP  )=FTC  15, NTP  HZC  J  HGR 
FT<  16»NTP  >=FT<  16, NTP  )+Z<  J  HRC  J  HGR 
B(  JJ  >=BC  JJ  >-ZX*<  COSA*BZ-SINA*DR ) 

GZ=  -ZX*C  COSA*BZ-SINA*DR  HTHETA/2.0 
FT<  5, NTP  )=FT<  5»NTP  HGZ 
FT<  6,  NTP  )=FTC  6, NTP  HRC  J  HGZ 
FT(  7»NTP  )=FTC  7 »NTP  )+Z<  J  HGZ 
FTC  8, NTP  )=FTC 8,  NTP  HZC  J  HRC  J  HGZ 
FTC  17»NTP  )=FT<  17»NTP  )+GZ 
FT<  18, NTP  )=FT(  18, NTP  HRC  J  HGZ 
FT<  19»NTP  >=FTC  19, NTP  HZC  J  HGZ 
FT(  20»NTP  )=FT<  20, NTP  HZC  J  HRC  J  HGZ 
640  CONTINUE 

701  IF<  NUMST .EG.O  )  GO  TO  700 
DO  680  L= 1» NUMST 
I=IT<  L  ) 

J=JTCL> 

RT=STCL>/6. 

RX=2,*RC  I  HRC  J  ) 

ZX=RC  I  H2.*RC  J) 

XX=SQRT<  <  RC  J  )-R<  I  )  H*2+C  ZC  J  )-Z<  I  )  H*2  ) 

II=3*I-KSHIFT 

JJ=3*J-KSHIFT 

IFC II.LE.O.OR.II.GT.ND)  GO  TO  670 
B< II  >=BC II  HRT*RX*XX 
GT=RT*RX*XX*THETA/2 . 0 
FT< 9, NTP  )=FTC 9,  NTP  HGT 
FT( 10»NTP  )=FTC 10, NTP  HRC I  HGT 
FT(  11, NTP  )=FT(  1 1 , NTP  HZC  I  HGT 
FTC  12, NTP  )=FTC  12, NTP  HZC  I  HRC  I  HGT 
FTC  21 » NTP  )=FTC  21 »NTP  HGT 
FTC  22»NTP  >=FTC  22, NTP  )+RC  I  HGT 
FTC  23, NTP  )=FTC  23, NTP  )+ZC  I  HGT 
FTC  24, NTP  >=FTC  24, NTP  )+ZC  I  HRC  I  HGT 
670  IFC JJ.LE.O.OR.JJ.GT.ND)  GO  TO  680 
BC  JJ  )=BC  JJ  HRT*ZX#XX 
GT=RT*ZX*XX*THETA/2 .0 
FTC  9 , NTP  )=FTC 9 , NTP  )+GT 
FTC  10, NTP  )=FTC  10, NTP  HRC  J  HGT 
FTC  1 1 , NTP  )=FTC  1 1 , NTP  )+ZC  J  )*GT 
FTC  12, NTP  >=FTC  12, NTP  HZC  J  HRC  J  HGT 
FTC  21 , NTP  )=FTC  21 ,  NTP  HGT 
FTC  22, NTP  )=FTC 22, NTP  HRC  J  HGT 
FTC  23, NTP  >=FTC  23, NTP  HZC  J  HGT 
FTC  24, NTP  )=FTC  24 , NTP  >+ZC  J  HRC  J  HGT 
680  CONTINUE 

c*.  a***********************  *********** 
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C  ADD  DISPLACEMENT  BOUNDARY  CONDITIONS 

C*  *********************  ************** 

'  700  DO  750  H=NL,NH 
I  DM=0 

IF! M.GT .NUMNP  )  GO  TO  750 
IF!  CODE( M  )»GT .3  )  GO  TO  751 

U=XR<  M  )  • 

N=3*M-2-KSHIFT 
752  IF(  CODE(  M  )  )  740,750,710 
710  I F(  CODE(  M  ) ♦ EQ » 1  )  GO  TO  720 
IF  !  CODE<  M  KEG. 2  )  GO  TO  740 
IF! CODEC  M  >.EG.3  )  GO  TO  730 
GO  TO  740 

720  CALL  MODIFY<  ND2»N,U  ) 

COD£< M  )=CODEC  M  )+IDM 
GO  TO  750 

730  CALL  MODIFY<  NB2»N»U  ) 

740  U=XZ!M> 

N=N+1 

CALL  MODIFY<  ND2,N»U  ) 

CODE<  M  >=CODE!  M  )+IDM 
GO  TO  750 
751  IDM=IDM+4 
U=XT( M  ) 

N=3*M-KSHIFT 

CALL  MODIFY! ND2*N,U) 

IJ=XR<  M  ) 

N=3*M-2~KSHIFT 
IF!  CODE! M  ) . EQ . 4  )  GO  TO  750 
CODE!  M  )=CCDE!  M  )-4 
GO  TO  752 
750  CONTINUE 

C*  *********************************** 
C  WRITE  BLOCK  OF  EQUATIONS  ON  FORTRAN  UNIT  AND  SHIFT  UP  LOWER  BLOCK 
C*  *********************************** 
WRITE  !  2  )  !B!N)»!A!N,M),M=1, MBAND  ) ,  N= 1 » ND  ) 

DO  800  N=1 »ND 

K=N+ND 

B!  N  )=B!  K  ) 

B!  K  >=0.00 
DO  800  M=1 »ND 
A!  N»M  )=A!  K,M  > 

800  A!  K ,  M  )=0 . 00 

C*  *********************************** 
C  CHECK  FOR  LAST  BLOCK 

C*  *********************************** 
IF! NM.LT .NUMNP  )  GO  TO  200 
IF! STOP.NE.O. >  STOP 

2000  FORMAT  ! 27H  NEGATIVE  AREA  ELEMENT  NO., 14) 

2001  FORMAT  ! 46H  BAND  WIDTH  EXCEEDS  ALLOWABLE  FOR  ELEMENT  NO., 14) 
RETURN 

END 

SUBROUTINE  STORE 
INTEGER  CODE 

C0MM0N/ANS4/  FT! 24,4  >,GTS1U! 24  >,GTS1UT!  24,4  > 

COMMON/NXDATA/NTP , NTS , NTOTS , GTS1G!  24 , 24 , 8 ) 

COMMON/NXMESH/THETANC  8  )»NPC!  8,8  ) 

COMMON/NPDATA/R!  10  ), CODE!  10  >,XR!  10  ),2!  10  >,XZ!10  ), 

1NPNUM!  4,  8  )»T! 10  >»XT! 10  ) 

COMMON/SOLVE/B!  72), A!  72,36 >,NUMBLK, MBAND 
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ounnUN/  ULDbtiU/ 1  it  i!4 » iJ  )  >  f  C(  24  i  y  ) » UC(  24 »  8  )  ,  SfU  24  » 24  » L  ) 
COMMON/BASIC/ ACELZ  *  ANGVEL,  ANGACC,  TREE  »  VOL*NUHii2  ,  KUi'.EI 

1NUMST 

COMMON/ ANS2/LVK  24  >»R1(  24,24  >»SK1(  24,24  > > 24,24  ) 
DIMENSION  MU< 24  ) 


1  {iUsIr'L  ,  UUi'i  liU  > 


DO  50  J=1 » 24  . 

50  RK  I,.J  >  =  0.0 
NS  =  NTP 
DO  HO  KK  =  1,4 
NP1  =  NPC<  NS»KK ) 

NP2  =  NPC< NS»KK+4 ) 

DO  110  1=  1,3 

RM3KKK-1  )+I  ,1*4-3  >  =  1.0 

Rl(  3*<  KK-1  HI  ,1*4-2  )  -  R(NP1> 

Rl<  3*<  KK-1  )+I  ,1*4-1  )  =  Z(NP1) 

Rl< 3*< KK-1  )+I  ,1*4  )  =  R<NP1)  *  Z<  NP1  ) 

Rl<3*<  KK-1  >+1+12,1*4+9  )  =  1.0 
Rl<  3*<  KK-1  )+I  +  12, 1*4+10)  =  R(NP2) 

Rl< 3*< KK-1 >+1+12,1*4+11 >  =  ZINP2) 

Rl<  3*<  KK-1  )+I  +  12, 1*4  +  12  )  =  R(  NP2  )*  Z(  NP2  ) 

UC<  3*(  KK-1  >+I»NS  )  =  B<3*NPl-3+I> 

UC<  3*<  KK-1  )+I  +  12  ,NS  )  =  B(3*NP2-3+I) 

110  CONTINUE 

CALL  MINU< R1,24,D1,LU,MU ) 

WRITE(  25  )  ( ( Rl(  I, J),J=  1,24  >,1-1,24) 

DO  115  1=1,24 
FE< I, NS  )=0.0 
FI( I, NS  )=0.0 
DO  115  J=1 ,24 

FE< I ,NS  )  =  FE< I ,NS  )+  Rl( J, I  )*FT< J,NTP  ) 

FI(  I ,  NS  )  =  FKI,  NS)  +  Rl<  J,  I  )*GTS1UT(  J,NTP  ) 

SK1< I»J  )  =  0.00 
DO  115  K=1 ,24 

SK1( I , J  )  =  SK  1< I , J  )  +  Rl<  K, I  )  *  GTS1G(  K, J,NTP ) 

115  CONTINUE 

DO  120  1=1,24 

DO  120  J=1 ,24 

SK( I , J,NS  )  =0.0 
no  120  K=1 f  24 

SK<  I  ,J,NS  )  =  SK<  I  ,J,NS )  +  SK1<I,K)  *  R1(K»J> 

120  CONTINUE 
RETURN 
END 

SUBROUTINE  STRESS 
INTEGER  CODE 

COMMON/VISC/EPSDN( 12, 10,8  >,SIGVP( 6 ),DEPSR( 6, 10,8 ),DELTIM 
COMMON/ ANS4/FT< 24,4  ),GTS1U<  24  )  >GTS1UT<  24,4  ) 

COMMON/BAS  IC/ACELZ ,  ANGVEL ,  ANGACC ,  TREF ,  VOL  ,NUMNP ,  N'UKEL  ,NUMPC ,  NUMSC , 
1NUMST 

COMMON/MATP/RO( 6 >,E< 12,16,6  ),EE( 16  ),A0FTS(6) 

COMMON/NPDATA/R<  10  ),CODE<  10  ),XR(10  ),Z(10  ),XZ(10  ), 

1NPNUM?  4,  8  ),T(  10  ),XT<  10  ) 

COMMON/ELDATA/BETA(  10  ),EPR<  10  ),PR(  4  ),SH<  4  >»IX<3  ,5), 

1  IP<  4  ),  JP<  4  ),  IS(  4  ),  JS<  4  )» ALPHA(  10  ),IT(4  ),JT<4  ), 

2ST< 4  ) 

COMMON/ARG/RRR<  5  >,ZZZ(  5 ),RR(  4  ),ZZ(  4  >,S<  15,15  ),P(  15  )»TT<  6  ), 

1H(  6, 15  ),CRZ<  6,6  ),XI(  10  ),ANGLE<  4  >,5IG<  18),EF'S(  18  ),N 
COMMON /NONAX I /Sl(  30,30  )»P1< 30  ), THETA, BS1(  6,30 ) 

COMMON/NXBATA/NTP , NTS, NTOTS , GTS1G(  24,24,3  ) 


Rl<  K, I  )  *  GTS1G(  K , J , NTP ) 


SK1<  I,K)  *  R 1  (  K ,  J  ) 
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COMMON/NXMESH/THETAN<  8 ) , NPC<  8 , 8 ) 

C0MMG?'I/ARG1/S1GH  18  )»EPSH  18  )»DEPSP<  12  )»CEPSP<  6,6  ) 

COMMON/SOLUE/B<  72)»A<  72,36  ).NUHBLK»MBANIi 
COMMQN/CONVRG/IBONE 

COMHON/RATE/DKPR , SIGPR ,  BVR » EVR , PSRATE<  10,8 ) » NR ATE 
COMMON/PL ANE/NPP 

COMMON/RESULT /BS<  6 . 1 5  )  ,  B<  6 . 6  ) » C<  6 » 6  ) , AR » BB<  6 » 9  )  *  CNS<  6,6) 

DIMENSION  LM(  4  ),TP<  6  )»TR<  3,3),Q<3) 

DIMENSION  QQ( 3  ) 

COMMON/DUM1 /SI TEM<  3 , 30  )  *  GTS 1<  24 , 24  )  >  TS<  6 , 24 ) 

DIMENSION  Pll<  24  ) 

C  ************************************ 

C  INITIALIZE 

C*  *********************************** 
REWIND  3 
XKE=0. 

XPE=0. 

MPRINT=0 
ERROR=  »005 
ID0NE=1 
NRATE=0 

DO  200  N=1 » NUMEL 
IX<  N»5  )=IABS<  IX<  N,5  )  ) 

CALL  QUAD 

MTYPE=IABS< IX(  N»5  ) ) 

DO  100  1=1.4 
11=3*1 

JJ=3*IX<  N.I ) 

PH  1 1-2)  =  B<  J.J-2  ) 

Pl<  II-l  )  =  B(  JJ-1 ) 

Pl<  II  )  =  B<JJ  ) 

PKII+10)  =  B<  JJ-2  ) 

Pl( 11  +  11  )  =  B< JJ-1) 

PH  11  +  12)  =  B<  JJ  ) 

P<  II-2)=B(  JJ-2) 

P<  II-l  )=B(  JJ-1  ) 

100  P<II)  =B<  JJ  ) 

READ<  3  )( <  GTSK  I .  J  ).  J=1 . 24  )» 1=1 » 24  ) 

DO  115  1=1.24 
GTS1U<  I  )=0.0 
DO  115  J=1 .24 

115  GTS1U<I)  =  GTS1U(  I  )+  GTS1<  I » J  )*P1(  J  ) 

DO  116  1=1.24 

116  GTS1UT< I.NTP  )=GTS1UT< I »NTP  )  +  GTS1U( I  ) 

DO  110  1=1,3 

110  Q<  I  )=P<  1+12  > 

DO  120  1=1,3 
DO  120  J=1 , 3 

120  TR(I,J  )=S(  I+12.J+12) 

CALL  SYMINV( TR»3 ) 

DO  125  J=1 .3 
QQ(  J  )=0,00 
DO  125  K=1 » 12 

OCX  J  )=QQ(  J  >+S(  J+12.K  )*P<  K  ) 

125  CONTINUE 

DO  130  1=1,3 
P<  1  +  12  )=0*00 
DO  130  J=1 ,3 

130  P<  1  +  12  )=P(  1+12  )+TR<  I ,  J  )*<  Q<  J  )-QQ<  J  ) ) 

500  CONTINUE 
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RETURN 

END 

SUBROUTINE  SYMINU! A. NMAX > 

'  DIMENSION  A( NMAX .NMAX) 

-  DO  300  N=1 .NMAX 

D=A<  N.N  ) 

DO  100  J=1 .NMAX 
100  A<  N.J  )=-A<  N»  J  )/D 
DO  210  1=1 .NMAX 
IF<  N.EQ. I  )  GO  TO  210 
DO  200  J=1 .NMAX 

IF!N.NE.J)  A(I»J  )=A<  I.J  )+A!  I,N)*A!N.J) 

200  CONTINUE 
210  A<  I  >N  >=A<  I  »N  )/D 
300  A<  N.N  )=1 .00/D 
RETURN 
END 

SUBROUTINE  TEMP(R.Z.T) 

COMMON/SOLVE/X! 888  ) . Y<  888  ) . TEM< 888  ) » NUMTC . MBAND 
DIMENSION  SMALL!  20  )» ISM<  20  ) 

C*  *********************************** 
C  INITIALIZE 

C*  *********************************** 
J=1 

JMAX=16 

IF(NUMTC.LT.JMAX  )  JMAX=NUMTC 
DO  10  1=1 .JMAX 
SMALL!  I  )=0. 

10  ISM(  I  )=0 

C*  *********************************** 
C  FIND  THE  JMAX  CLOSEST  POINTS 

C*  *********************************** 
DO  50  1=1. NUMTC 
DSQ=(  X<  I  )-R  )**2+(  Y<  I  )-Z  )**2 
IF<  DSO.GT • • IE-4  )  GO  TO  20 
T=TEM<  I  ) 

RETURN 

20  IF(I.EQ.l)  SMALL!  1  )=DSQ 
IF!  I.EQ.l  )  ISM!  1  )=1 
IF!  I »EQ. 1 )  GO  TO  50 

IF!  SMALL!  J  ) . LE . DSQ . AND . J . LT . JMAX  )  SMALL! J+l )=DSQ 
IF!  SMALL!  JKLE.DSG.  AND.  J.LT.  JMAX)  ISM!  J+l  )=I 
IF! SMALL! J  ).LE.DSQ  )  GO  TO  40 
DO  30  K=1.J 
JB=J-K  +1 

IF!  JB.EQ.O  )  GO  TO  40 
SMALL!  JB+1  )=SMALL!  JB  ) 

ISM!  JB+1  )=ISM!  JB) 

SMALL! JB  )=DSQ 
ISM!  JB  )=I 

IF! JB.EO. 1 )  GO  TO  40 
IF!  SMALL!  JB-1  ).LE.  DSQ)  GO  TO  40 
30  CONTINUE 
40  IF!  J.LT. JMAX)  J=J+1 
50  CONTINUE 

C*  *********************************** 
C  FIND  THE  THIRD  TEMPERATURE  POINT  BY  AREA  TEST 

C*  *********************************** 
JCHK=JMAX-2 
J=0 


c 

C  MATRIX  P  NOW  CONTAINS  15  DISPLACEMENTS  FOR  QUADRILATERAL  ELEMENT 
C 

C'  CALCULATE  AVERAGE  STRAINS 

e 

DO  140  1=1,6 
EPS<  I  )=0  *00 
DO  140  J=1 *  15 

140  EPS*  I  )=EPS(  I  )+BS<  I ,  J  )*P*  J  > 

C 

C  CALCULATE  AVERAGE  STRESSES 

C 

DO  151  1=1 ,6 

S1G<  I  )=EPSDN<  I »N,NTP  > 

DO  151  J=l,6 

151  SIG(  I  >=SIG<  I  )+CRZ<  I » J  )*EPS<  J  ) 

DO  152  1=1,6 
152  SIG<  I  )=SIG(  I  )-TT<  I ) 

C 

C  CALCULATE  STRAINS  IN  N-S-T  COORDINATES 

C 

DO  150  1  =  1,6 
EPS*  1+6  )=0.00 
DO  150  J=l»6 
DO  150  K=1 ,6 

150  EPS*  1+6  )=EPS*  1+6  )+D<  1 » J  )*C*  J»K  )*EPS*  K  ) 

C 

C  CALCULATE  STRESSES  IN  N-S-T  COORDIATES 

C 

DO  160  1=1 ,6 

SIG* 1+6  )=EPSDN< I+6»N»NTP  ) 

DO  160  J=1 *6 

160  SIG*  1+6  )=SIG(  1+6  )+CNS<  I  ,  J  )*EPS<  J+6  ) 

DO  161  M=1 » 6 
P*M)=0.00 
DO  161  II=1» 3 

IF*  AOFTS* MTYPE  )»EQ.l .  )  P<  M  )=CNS<  M ,  T  T  )*EE* II+9  > 

161  P*M)=P*M)+*T*N)-TREF)*CNS*M»II  )*EE<  II+?) 

DO  162  1=1,6 

162  SIG*  1+6  )=SIG<  1+6  )-P*  I  ) 

C 

C 

DO  300  1=1,12 
300  EPS*  I  >=100.0*EPS*  I) 

IF*  MPRINT ,NE ,0  )  GO  TO  210 
WRITE*  6,2000) 

WRITE* 6,2002) 

MPRINT=19 

210  MPRINT=MPRINT-1 

WRITE*  6,2001 )  N,RRR(  5  ),ZZZ*  5  ),<  SIG*  I ),  1=1, 12  ) 

WRITE*  6,2003  )  T*  N  ),<  EPS*  I ),  1  =  1 , 12  ) 

200  CONTINUE 

2000  FORMAT*  129H1  EL  R  Z  SIGMAR  SIGMAZ  SIGMAC  SIGMA 

1RZ  SIGMAZC  SIGMACR  SIGMAN  SIGMAS  SIGMAT  SIGMAN3  SIGMAST 

2  SIGMATN  ) 

2001  FORMAT* 1H0, 15, 1X»2F7 , 4 » 12F? . 0  ) 

2002  FORMAT*  128H0  TEMPERATURE  EPSR  EPSZ  EPSC  EPSR 

1Z  EPSZC  EPSCR  EPSN  EPSS  EPST  EPSNS  EPSST 

2  EPSTN ) 

2003  FORMAT* 6X,F13.0,2X, 12F7, 5  ) 
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h=ism<  i )  . 

1 2=  ISM:  2) 

.  60  I3=ISM<J+3> 

AREA- .  50*!  Y!  II  >*X!  13  )-Y(  I3)*X!  II  )  LV!  13  >*X<  12  )-Y<  I2)*X(  13  >  + 

'  1  Y<  12  )*X<  II  )-Y<  II  )*X!  12) ) 

Dl=<  X!  12  )~X!  II  >  >**2i!  Y!  12  >-Y(  II  )  )**2 
C  IF  Dl  IS  APPROXIMATELY  0.  IT  IS  ASSUMED  THAT. THERE  EXISTS  A 
C  DUPLICATION  OF  INPUT 

IFCB1.GT..1E-3)  GO  TO  70 

12=13 

J=J+1 

GO  TO  60 

70  IF< AREA**2.GT. ,1*U1*SMALL! 1  ))  GO  TO  SO 
J=.J+1 

IF!  J.LT.JCHK)  GO  TO  60 
WRITE!  6f2000  >  I1fI2fI3fJ 
T=TEM!  II  ) 

RETURN 

C*  ******************  *  *  *  ************** 
C  FIND  TEMPERATURE  INTERCEPT 

C************************************ 
80  DETA=Y<  II  >*!  TEM(  13  )-TEM<  12  >  >+Y!  12  )*!  TEM(  II  )-TEM!  13  )  ) 

1  +Y<  13  >*!  TEM(  12  )-TEM<  II)) 

DETB=X<  II  )*<  TEM(  I2)-TEM<  I3))+X(  12)*!  TEM(  I3)-TEM(I1  )) 

1  +X(  I3)*<TEM<  II  >-TEM<  12  >) 

DETC=TEM!  II  )*(  X(  12  >*Y<  13  )-X(  13  )*Y(  12  )  >+TEM<  I2)*(  X(  13  )*Y<  II  >-X(  ID* 
1  Y(  13  )  )+TEM(  13  )*<  X<  II  )*Y<  12  )-X<  12  )*Y<  II ) ) 

T=< DETA*R+DETB*Z+DETC  )/( 2.* ARE A  ) 

2000  FORMAT  < 28H  ERROR  IN  TEMPERATURE  INPUTfSH  11=14 »5H  12=14 f 

15H  I3=I4» 4H  J=I4 ) 

RETURN 

END 

SUBROUTINE  TEM2< NUMNP ) 

INTEGER  CODE 

COMMQN/NPDATA/R<  10  )rCOD£<  10  )»XR(  10  )»Z<  10  )»XZ(10  )f 

1NPNUM!  4 »  8  )fT!  10  )»XTC10  ) 

READ!  5  f  1000)  TCONST 
DO  100  N= If NUMNP 
100  T! N  )=TCONST 
1000  FORMAT! F10.0  X 
RETURN 
END 

SUBROUTINE  TRISTF  !IIfJJfKK> 

INTEGER  CODE 

COMMON/VISC/EPSDN!  12 f 10  f  8  )» SIGVP!  6  JfDEPSR! 6f10f8  >f DELTIM 
C0MM0N/PLAS/ALFA!6f  4  f  8  )  f  S I GYLD( 7  f  6  f  8  )  f IFGPL<  4f8) 
COMMON/NXDATA/NTPfNTSfNTOTSfGTSIG! 24f24»8  ) 

C0MM0N/MATP/R0(6)fE(  12f16f6)fEE(  16)fA0FTS(6) 

COMMON/BASIC/ACELZ  f  ANGVEL  f  ANGACC  f  TREF  f  VOLfNUHNP  f  NUMEL  f  NUMPC  f  NUMSC  f 
1NUMST 

C0MM0N/ARG/RRR(5)fZZZ<5)fRR(4)fZZ(4)fS(  15»15)fF'U5  >fTT<  6)f 
1H<  6f  IS  )fCRZ(  6f6  )fXI<  10  )f  ANGLE!  4  )fSIG!  18)»EPS(  18  >fN 
COMMON/NPDATA/R!  10  JfCODE!  10  )fXR(  10  )fZ(10  )fXZ<10  )f 

1NPNUM!  4 f  8  )»T<  10  )fXT!  10  ) 

COMMON/ELDATA/BETA!  10  )fEPR!  10  )fPR<4  )fSH(  4  )fIX!8  f5)f 
1  IP!  4  )f  JP!  4  )f  IS!  4  )fJS!4  )fALPHA(10  >fIT<4  )fJT(  4  >f 
2ST! 4  ) 

C0MM0N/N0NAXI/S1!  30  f  30  )  f  PI! 30  )  f  THETAf  BS1! it  30 ) 

COMMON/RESULT/BS!  6f15  )fD<  6f6  )fC!  6»6  >fARfBB!6f9  )fCNS(  6f6  ) 

DIMENSION  B1A!6f9  JfBIB!  6 f  9  )fB2A!  6 f 9  >fB2B!  6  f 9  >  fB3A!  6  r 9  )  f  B3B!  6f9) 
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DIMENSION  Bl<  6,9  )»B2<  6,9  ),B3<  6,9  >*F<  6,9  )»G<  9,6  >»V<  9;9  > 

DIMENSION  BF<  3  > , BFR(  3  ) » BFZ<  3  )  •  TF‘(  9  ) ,  B(  9 » 9  )  *  TVF'(  9  ) 

MTYPE=IABS< IX(  N»5  ) ) 

RR<  1  )=RRR(  II) 

RR<  2  )=RRR<  JJ  ) 

RR<3)=RRR(NN) 

ZZ<  1  >=ZZZ<  II  ) 

ZZ(  2  )=Z2Z(  JJ  ) 

ZZ<  3  )=ZZZ<  KK  > 

CALL  INTER 
VOL=VOL+XI< 1  ) 

COMM=RR<  2  >*<  ZZ<  3  >-ZZ<  1 )  )+RR(  1  >*<  ZZ<  2  )-ZZ<  3 )  )+RR(  3  >*(  ZZ<  1  )-ZZ(  2  )  ) 
DO  10  1  =  1,6 
DO  10  J=1 ,9 
Bl<  I ,  J  >=0.00 
B2<  I » J  )=0»00 
B3<  I » J  )=0.00 

FILL  B1  MATRIX-CONSTANT  TERMS 
B1C  1 , 1  )=(  ZZ<  2  >-ZZ<  3  >  )/COMM 
Bl(  1,4  )=(  ZZ<  3  >-ZZ<  1  )  )/COMH 
Bl(  If  7  )=<  ZZ(  1  )-ZZ<  2  )  )/COMM 
Bl<  3,1  )=B1(  1,1) 

Bl<  3,4  )=B1<  1 » 4  ) 

Bit  3,7  )=B1<  1 » 7  ) 

Bl(  2,2  )=<  RR(  3  )-RR(  2  )  )/COMM 
Bl<  2,5  )=(  RR(  1  )-RR(  3  )  )/COMM 
Bl(  2,8  )=<  RR(  2  )-RR(  1  )  >/COMM 
B1<4»1  )=B1(  2»  2  ) 

Bl(  4,4  )=B1(  2,5  ) 

Bl<  4,7  )=B1<  2,8 ) 

Bl(  4,2  )=B1<  1 » 1  ) 

B1(4»5)=B1<  1,4) 

Bl<  4,8  )=B1<  1,7) 

Bl<  5,3  )=B1(  4,1) 

Bl<  5,6  )=B1<  4,4 ) 

Bl<  5,9  )=B1<  4,7  ) 

FILL  B2  MATRIX-l/R  TERMS 

B2<  3,1  )=<  1/COMM  >*<  <  ZZ<  3  >-ZZ<  2  )  )*RR(  2  >+(  RR(  2  )-RR(  3  )  >*ZZ<  2 ) ) 

B2<  3,4  )=<  1/COMM  )*<  <  ZZ<  1  )-ZZ(  3  >  >*RR<  3  )-(  RR<  1  )-RR<  3  )  )*ZZ(3 )  ) 

B2<  3,7  )=<  1/COMM  >*<  <  ZZ(  2  )~ZZ(  1  )  >*RR<  1  >+<  RR(  1  )-RR<  2  )  >*ZZ<  1 )  ) 
B2(6,3)=-B2<3»1  ) 

B2<  6,6  )=-B2(  3,4  ) 

B2(  6 » 9  )=-B2(  3 » 7  ) 

FILL  B3  MATRIX-Z/R  TERMS 
B3<  3, 1  )=<  RR<  3  )-RR<  2  )  >/COMM 
B3<  3,4  >=(  RR<  1  )-RR(  3  )  J/COMM 
B3<  3,7  )=<  RR<  2  )-RR<  1  )  )/COMM 
B3(  6,3  )=(  RR<  2  )-RR(  3  )  )/COMM 
B3<  6,  6  )=(  RR(  3  )-RR<  1  )  )/COMM 
B3(  6,9  )=<  RR<  1  )-RR(  2  )  )/COMM 
AR=AR+XI( 1  ) 

DO  80  1=1*6 
DO  80  J=l,9 

BB<  I  *  J  )=B1<  I ,  J  >*XI<  1  )+B2<  I » J  >*XI<  2  >+B3<  I » J  )*XI<  4  ) 

DO  81  K=1 ,6 
DO  81  1  =  1,3 

BS(  K,3#JJ-3+I  )=BB< K,  1+3  >+BS<  K,3*JJ-3+I  ) 

BS<  Kt 3*11-3+1  >=BB<  K»  I  >+BS(K, 3*11-3+1 ) 

BS<K,3*KK-3+I  )=BB( K, 1+6  )+BS<  K»3*KK-3+I  ) 

DO  220  1=1,6 


UU  *20  J=l»9 

B1A(  I » J  )=B1(  I,J  »fc>:i<  1  >+B2(  I,J  )*XI<  2  )  i B3<  I , J  >*XI<  4  ) 

B2A<  I*  J  >=B1<  I,  J  )*XI(  2  )+B2(  It  J  >*XI<  3  >+B3<  I*  J  >*XI<  5  ) 

B3A(  I»J  )=B1<  I»J  )*XI<  4)+B2(  I ,  J  >*XI<  5  >+B3<  I,  J  )*XI<6> 

*  220  CONTINUE 

DO  200  1  =  1,6 
BO  200  K=l»9 
B1B(  I  »K  )=0.0 
B2B(  I  f  K  )=0*0 
B3B<  I,K>=0.0 
DO  200  J=1 , 6 

B1B(  I? K  )=B1B(  I f K  >+CRZ<  I»J  >*B1A<  J,K> 

B2B<  I  f  K  )=B2B(  I  f  K  )+CRZ<  I  f  J  )*B2A<  J  >  K  ) 

B3B(  IfK)=B3B(  IfN  )+CRZ<  I,J  )*B3A<  J,K) 

200  CONTINUE 

DO  230  1=1 f 9 
DO  230  K= 1 f  9 
B(  If K  >=0.0 
DO  230  J=1 f  6 

B<  I  »K  )=B<  I  f  K  )+Bl<  J » I  )*B1B<  J»K)+B2<  J,I  >*B2B<  J,K)+B3<  J,I  >*B3B<  JfK) 
230  CONTINUE 

C  ASSEMBLE  QUADRILATERAL  STIFFNESS  MATRIXf  Sf  FROM  TRIANGULAR 
C  STIFFNESS  MATRIXf  B. 

IIM=3*II-3 
JJM=3*JJ-3 
KKM=3*KK-3 
DO  120  1=1 f 3 
DO  120  J=1 f 3 

S< IIM+I f IIM+J  )=B< I  fj  >+S< IIM+I, IIM+J > 

S<  IIM+I fJJM+J  )=B<  I  f  J+3  )+S<  IIM+I,  JJM+J) 

S< IIM+I f KKM+J  )=B< I  , J+6 >+S( IIM+I, KKM+J) 

S(  JJM+I, IIM+J  )=B(  1+3,  J  )+S( JJM+I, IIM+J  ) 

S(  JJM+I fJJM+J  >=B( 1+3 , J+3  )+S( JJM+I , JJM+J  ) 

S<  JJM+I , KKM+J  )=B( 1+3, J+6 >+S<  JJM+I , KKM+J ) 

S<  KKM+I ,  IIM+J  )=B<  1+6,  J  )+S(KKM+I,  IIM+J) 

S<  KKM+I fJJM+J  )=B< 1+6 , J+3  )+S<  KKM+I, JJM+J  ) 

S<  KKM+I , KKM+J  )=B( 1+6, J+6  )+S<  KKM+I , KKM+J ) 

120  CONTINUE 

C  ASSEMBLE  BODY  FORCES  MATRIX 

BF(  1  )=(  ZZ(  3  )*RR<  2  )-RR<  3  )*ZZ<  2  )  >/COMM 
BF<  2  )=<  ZZ(  1  )*RR<  3  )-RR(  1  >*ZZ<  3  )  )/COMM 
BF<  3  )=<  ZZ(  2  )*RR<  1  >-RR<  2  )*ZZ<  1  )  >/COMM 
BFR<  1  )=(  ZZ<  2  )-ZZ<  3  )  )/COMM 
BFR<  2  )=(  ZZ(  3  )-ZZ<  1  )  )/COMM 
BFR<  3  )=<  ZZ<  1  >-ZZ<  2  )  )/COMM 
BFZ(  1  )=(  RR(  3  )-RR<  2  )  )/COMM 
BFZ<  2  )=<  RR<  1  )-RR(  3  )  )/COMM 
BFZ(  3  )=<  RR<  2  )-RR<  1  )  )/COMM 
C  BODY  FORCE  IN  Z-DIRECTION 
COMM=-ACELZ*RQ< MTYPE  ) 

DO  140  1=1,3 
IIK=3*I-1 

140  TP<  IIK  )=COMM*<  BF(  I  )*XI<  1  )+BFR<  I  )*XI<  7  )+BFZ<  I  )*XI<  8  ) ) 

C  BODY  FORCE  IN  R-DIRECTION 
C0MM=ANGVEL**2*R0<  MTYPE  ) 

DO  150  1=1,3 
L=3*I-2 

150  TP<  L  )=COMM*<  BF<  I  )*XI<  7  )+BFR<  I  )*XI<  9  )+BFZ<  I  >*XI<  10  > ) 

C  BODY  FORCES  IN  TANG.  DIRECTION 
COMM=-ANGACC*RO<  MTYPE  > 
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IIM=3*I 

160  TP<  IIM  )=C0MM*<  BF<  I  >*XI<  7  >+BFR<  I  >*XI'.  9  >+BFZ<  I  >*XI<  10  ) ) 

ADD  THERMAL  EFFECTS 
DO  161  J=l»9 
DO  161  K=1 »6 

161  TP(  J)=TP<  J)+<XI<l)*Bl<K,J)+XI<2)*B2<KyJ) 

11X1(4  )*B3<  K»  J  )  )¥TT(  K  ) 

REARRANGE  TP  INTO  P-MATRIX. THE  BODY  FORCES  MATRIX 

K=3*I 1-2 

L=3*JJ-2 

M=3*KK-2 

DO  170  1  =  1,3 

J=T-1 

Pl(KtJ  )  =  P1<K+J  >+TP<  I  >*THETA/2.0  +TP<  1+6 >*THETA/4 . 0 
PKK+J+15)  =  Pit  K+J  +  15  )+TP<  I  >*THETA/2.0+TP<  1+6  >*THETA/4 .0 
P1<L+J  )  =  PHL+.J  )+TP<  1+3  )¥THETA/2*0+TP<  1+6  )¥THFTA/4.0 
Pl<  L+J+15  >  =  Pl(  L+J+15  )+TF‘<  1+3  )*THETA/2. 0+TP<  1+6  )¥THETA/4»0 
P1<M+J  )  =  PKM-hJ  )+TP(  1+6  )¥THETA/2*0 
PKM+J+15)  =  Pl<  H+J+15  )+TP<  1+6  )*THETA/2.0 
P<  K+J  >=P<  K+J  )+TP<  I  ) 

P<L+J)=P<L+J)+TP<  1+3) 

170  P<  M+J  )=P<  M+J  >+TP<  1+6  > 

IF(  IFGPL<  N»NTP  ).EQ.O  )  GO  TO  190 
DO  174  1=1 » 9 
TVP<  I  )=0.0 
DO  174  J=l»6 

TVP<  I  )=TVP<  I  )+BB<  J » I  >*EPSDN<  J, N»NTP  ) 

174  CONTINUE 

DO  180  1=1 »3 
J=I-1 

P<  K+J  )=P<  K+J  )-TVP<  I  ) 

P(  L+J  )=P<  L+J  )-TVP(  1+3  ) 

180  P<  M+J  )=P(  M+J  )-TVP<  1+6  ) 

190  CONTINUE 
RETURN 
END 

SUBROUTINE  XMODFY<U,N> 

COMMON/NXSOLV/SK  <36  ,24  >,R1  < 132 >,FTOT< 132 >,NSZF 
NBAND=24 
DO  10  M=2,NBAND 
K=N-M+1 

IF<  K *LE»0  )  GO  TO  5 
R1<K)=R1(K)-SN<K,M)*U 
SK(  K»M  >=0» 

5  K=N+M-1 

IF<  NSZF.LT .K >  GO  TO  10 
Rl<  K  )=R1<  K  )-SK(  N»M  )*U 
SK<  N»M  )=0. 

10  CONTINUE 
SK'<  N » 1  )=  1  * 

Rl<  N  )=U 

RETURN 

END 

SUBROUTINE  XSOLVE 

COMMON/NXSOLV/SK  <36  ,24)>R1  <  132  ),FTOT<  132  ),NSZF 
NBAND=24 
DO  300  N=1 »NSZF 
I=N 

DO  290  L=2»NBAND 


1=1+1 

IF<SK(N,L))  240 *290 >240 
,  240  AC=SK<  N,L  )/SK(  {-4*1) 

'  J=0 

"  DO  270  i\=L > NBAND 

J=J  +  1 

IF<  SK( N ,K  ) )  260 1 270 » 260 
260  SK<  1 1 0  )=SK(  I » -J  )-AC*SK(  N ,  K  ) 

270  CONTINUE 
280  SK<  N»L  )=AC 
C 

Rl<  I  )=R1(  I  )-AC#Ri(  N  > 

290  CONTINUE 

300  Rl<  N  )=R1<  N  )/SK<  N>  1  ) 

C 

N=NSZF 
350  N=N-1 

IF<  N  )  500,500,360 
360  L=N 

DO  400  K=2» NBAND 
L=L+1 

IF<  SN(  N»K  ) )  370,400,370 
370  R1<N)=R1(N)-SK<N,N)*R1<L) 

400  CONTINUE 
GO  TO  350 
C 

500  RETURN 
ENO 

SUBROUTINE  YIELB(  N ,  NS , MTYPE  ) 

DIMENSION  DALFA<  6  )  » SIGYB<  3  ) 

COMMQN/PLAS/ALFA<  6,  4,8  )*SIGYLB(  7,6,8  )»IFGFL(  4,8) 

COMMON/ INCR/NOL INC , NOL , I NERT » NUMMAT  »  S I GTOT(  1 2 ,  4,8) 

1 ,EPSTOT< 12,  4,8) 

COMMON/ ARGl/SIGK  18),EPS1(  18),DEPSP<  12  ),CEPSP(  6,6  ) 

C=SIGYLD<  7 » MTYPE , NS ) 

DO  50  1=1,6 

50  ALFA(  I , N , NS  )=ALFA< I , N , NS  )+C*DEPSP<  I +6 ) 

C  WRITE( 6,1000 )N,NS 

C1000  FORMAT< "  ALFA  FOR  EL  "»I5,"  SEGMENT",  15) 

C  WRITE<  6,1100  )<  ALFA(  I  ,N,NS  ),  1  =  1 ,6  ) 

1100  FORMATC  ",6E12.6) 

1  DO  100  1=1,6 

100  5IGK  I  )=SIGTOT<  I  V6 ,N, NS  )-ALFA(  I,N.NS  ) 

C  GET  COMBINATION  YIELD  STRESSES 

SIGYB< 1  )=1 ./SIGYLDC 1, MTYPE, NS  )**2-l . /SIGYLD< 2, MTYPE, NS )**2 
1  -l./SIGYLD<3,MTYPE,NS)**2 

SIGYB< 2  )=  1 ,/SIGYLIK  2, MTYPE, NS  )**2-l . /SIGYLEK  1 , MTYPE, NS )**2 
1  -1,/SIGYLD< 3,MTYPE,NS)**2 

SIGYB(  3  )  =  1  ,/SIGYLD<  3, MTYPE, NS  5**2- 1  ./SIGYLEK  2, MTYPE, NS  )U2 
1  -l./SIGYLDt 1, MTYPE, NS  )**2 

C  TEST  YIELD  CRITERION  ********* 

TEST=0.0 
DO  200  1=1,6 

200  TEST=TEST+SIG1(  I  )#*2/SIGYLD<  I, MTYPE, NS  >U2 

TEST=TEST+  SIGYB<  1  >*SIG1<  2  )*SIG1(  3  )  +SIGYEK  2  )*3IG1<  1  )*SIG1(  3  ) 
1  +SIGYB<  3>*SIG1<  1  )*SIG1<  2) 

IFGPL<N,NS  )=0 

IF  <  TEST *LE« 1 ,0  )  GO  TO  500 
IFGPUN,NS  )=1 
500  RETURN 
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